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ABSTRACT 


This  work  accomplishes  signal  denoising  using  the  Bootstrap  method  when  the 
additive  noise  is  Gaussian.  The  noisy  signal  is  separated  into  frequency  bands  using  the 
Fourier  or  Wavelet  transform.  Each  frequency  band  is  tested  for  Gaussianity  by 
evaluating  the  kurtosis.  The  Bootstrap  method  is  used  to  increase  the  reliability  of  the 
kurtosis  estimate.  Noise  effects  are  minimized  using  a  hard  or  soft  thresholding  scheme 
on  the  frequency  bands  that  were  estimated  to  be  Gaussian.  The  recovered  signal  is 
obtained  by  applying  the  appropriate  inverse  transform  to  the  modified  frequency  bands 
The  denoising  scheme  is  tested  using  three  test  signals.  Results  show  that  EFT- based 
denoising  schemes  perform  better  than  WT- based  denoising  schemes  on  the  stationary 
sinusoidal  signals,  whereas  WT-based  schemes  outperform  FFT-based  schemes  on 
chirp  type  signals.  Results  also  show  that  hard  thresholding  never  outperforms  soft 
thresholding;  at  best  its  performance  is  similar  to  soft  thresholding. 


v 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


vi 


TABLE  OF  CONTENTS 


I.  INTRODUCTION . 1 

II.  BACKGROUND . 3 

A.  PROCESSING  TECHNIQUES . 3 

1 .  F ourier  Analysis . 3 

a  Fourier  Series . 3 

b.  Fourier  Transform . 4 

c.  Short-Time  Fourier  Transform . 5 

2.  W  avelet  Analysis . 6 

a.  The  Continuous  Time  Wavelet  Transform . 7 

b.  The  Discrete  Time  Wavelet  Transform . 7 

c.  Mallat’s  Algorithm . 8 

B.  HIGHER  ORDER  STATISTICS . 8 

1.  Moments  and  Cumulants . 8 

2.  Variance,  Skewness  and  Kurtosis  Measures . 9 

C.  THE  BOOTSTRAP . 10 

1.  Basic  Principle . 10 

2.  Parameter  Confidence  Interval . 11 

3.  Hypothesis  Testing . 11 

a.  Description . 11 

b.  Estimation  of  the  Standard  Deviation  for  p . 13 

III.  DENOISING  USING  THE  KURTOSIS  AND  THE  BOOTSTRAP . 15 

A.  INTRODUCTION . 15 

B.  FAST  FOURIER  TRANSFORM  BASED  DENOISING . 15 

1.  Short-Time  Fourier  Transform  Step . 16 

2.  Gaussianity  Test  Step . 16 

3.  Thresholding  Step . 16 

a.  Hard  Thresholding . 17 

b.  Soft  Thresholding . 17 

4.  Inverse  F ourier  T ransform . 18 

C.  WAVELET  TRANSFORM  BASED  DENOISING . 19 

1.  The  Wavelet  Transform . 19 

2.  Gaussianity  Test . 19 

3.  Thresholding . 20 

4.  Inverse  W avelet  T ransform . 20 

IV.  SIGNAL  DESCRIPTION  AND  SIMULATION  RESULTS . 21 

A.  SIGNAL  DESCRIPTION . 21 

1.  Sinusoidal  Signal . 21 

2.  Constant  Amplitude  Chirp . 21 

3.  Chirp  With  an  RC  Time  Constant -Like  Amplitude  Increase . 22 

vii 


B.  SIMULATION  RESULTS . 22 

1 .  Fast  F ourier  T ransform-Based  Denoising . 24 

a.  Sinusoidal  Signal . 24 

b.  Constant  Amplitude  Chirp . 25 

c.  Chirp  With  an  RC  Time  Constant-Like  Amplitude 

Increase . 26 

2.  Wavelet  Transform-Based  Denoising . 27 

a.  Sinusoidal  Signal . 28 

b.  Constant  Amplitude  Chirp . 29 

c.  Chirp  With  an  RC  Time  Constant-Like  Amplitude 

Increase . 30 

V.  CONCLUSIONS  AND  FUTURE  WORK . 33 

A.  CONCLUSIONS . 33 

B.  FUTURE  WORK . 34 

APPENDIX  A.  MATLAB  CODES . 35 

A.  FFT-BASED  DENOISING . 35 

B.  WT-B  ASED  DENOISING . 37 

C.  SUB  FUNCTIONS . 49 

1.  MakeSignal . 49 

2.  Kurtosistest . 55 

APPENDIX  B.  SIMULATION  RESULTS . 57 

A.  FFT-BASED  DENOISING . 57 

1.  Sinusoidal  Signal . 57 

2.  Constant  Amplitude  Chirp  Signal . 62 

3.  Chirp  with  an  Amplitude  that  Increases  in  an  RC  Time 

Constant  Fashion . 66 

B.  WT-B  ASED  DENOISING . 70 

1.  Sinusoidal  Signal . 70 

2.  Constant  Amplitude  Chirp . 75 

3.  Chirp  with  an  Amplitude  that  Increases  in  an  RC  Time 

Constant  Fashion . 79 

APPENDIX  C.  COMPARISON  WITH  WAVELET  SHRINKAGE 

ALGORITHM . 83 

A.  WAVELET  SHRINKAGE  ALGORITHM . 83 

1.  Choosing  a  Threshold  Value . 83 

a.  Stein’s  Unbiased  Risk  Estimator  (SURE)  Threshold. . 83 

b.  Sqtwolog  Threshold . 83 

c.  Heursure  Threshold. . 83 

d.  Minimax  Threshold . 84 

2.  Thresholding  Methods . 84 

a.  Hard  Thresholding . 84 

b.  Soft  Thresholding . 84 

B.  COMPARISONS . 84 

1 .  Sinusoidal  Signal . 85 

viii 


2.  Constant  Amplitude  Chirp . 86 

3.  Chirp  With  an  RC  Time  Constant -Like  Amplitude  Increase . 87 

LIST  OF  REFERENCES . 89 

INITIAL  DISTRIBUTION  LIST . 91 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


x 


LIST  OF  FIGURES 


Figure  1.  (a)  Time -Frequency  resolution  for  the  STFT,  (b)  Time -Frequency 

resolution  for  the  WT . 6 

Figure  2.  Mallat’s  algorithm . 8 

Figure  3.  FFT-based  denoising  scheme . 15 

Figure  4.  WT-based  denoising  scheme . 19 

Figure  5.  Sinusoidal  test  signal;  frequency  /0  =  1 0/5  1 2 ,  sampling  frequency  fs  =  1...21 

Figure  6.  (a)  Second  test  signal,  (b)  Third  test  signal . 22 

Figure  7.  Distance  measure  d  and  MSE  performance  criteria  for  the  sinusoidal 

signal;  FFT-based  denoising . 25 

Figure  8.  Distance  measure  d  and  MSE  performance  criteria  for  constant  amplitude 

chirp  signal;  FFT-based  denoising . 26 

Figure  9.  Distance  measure  d  and  MSE  performance  criteria  for  increasing 

amplitude  chirp  signal;  FFT-based  denoising . 27 

Figure  10.  Distance  measure  and  MSE  for  test  signal  1,  using  a  Wavelet-based  4- 
level  decomposition  and  method  2,  with  Daubechies  wavelet  orders  3,  4 

and  5 . 29 

Figure  11.  Performance  for  constant  amplitude  chirp  signal,  using  a  Wavelet  based  4 
level  decomposition  and  method  1,  with  Daubechies  wavelet  orders  3,  4 

and  5 . 30 

Figure  12.  Performance  for  increasing  amplitude  chirp  signal,  using  a  Wavelet-based 
4-level  decomposition  and  method  1,  with  Daubechies  wavelet  orders  3,  4 

and  5 . 31 

Figure  13.  FFT-based  denoising  for  sinusoidal  signal  type . 58 

Figure  14.  FFT-based  denoising  for  sinusoidal  signal  type,  SNR=-6  dB . 59 

Figure  15.  FFT-based  denoising  for  sinusoidal  signal  type,  SNR=0  dB . 60 

Figure  16.  FFT-based  denoising  for  sinusoidal  signal  type,  SNR=6  dB . 61 

Figure  17.  FFT-based  denoising  for  constant  amplitude  chirp  type  signal . 62 

Figure  18.  FFT-based  denoising  for  constant  amplitude  chirp  type,  SNR=-6  dB . 63 

Figure  19.  FFT-based  denoising  for  constant  amplitude  chirp  type,  SNR=0  dB . 64 

Figure  20.  FFT-based  denoising  for  constant  amplitude  chirp  type,  SNR=6  dB . 65 

Figure  21.  FFT-based  denoising  for  increasing  amplitude  chirp  type  signal . 66 

Figure  22.  FFT-based  denoising  for  increasing  amplitude  chirp  type,  SNR=-6  dB . 67 

Figure  23.  FFT-based  denoising  for  increasing  amplitude  chirp  type,  SNR=0  dB . 68 

Figure  24.  FFT-based  denoising  for  increasing  amplitude  chirp  type,  SNR=6  dB . 69 

Figure  25.  WT-based  denoising  for  sinusoidal  signal  type,  4  level,  method  2,  with 

Daubechies  wavelet  orders  3,  4  and  5 . 71 

Figure  26.  WT-based  denoising  for  sinusoidal  signal  type,  4  level,  method  2, 

Daubechies  wavelet  4,  SNR=-6  dB . 72 

Figure  27.  WT-based  denoising  for  sinusoidal  signal  type,  4  level,  method  2, 

Daubechies  wavelet  4,  SNR=0  dB . 73 


xi 


Figure  28.  WT-based  denoising  for  sinusoidal  signal  type,  4  level,  method  2, 

Daubechies  wavelet  4,  SNR=6  dB . 74 

Figure  29.  WT-based  denoising  for  constant  amplitude  chirp  type,  4  level,  method  1, 

with  Daubechies  wavelet  orders  3,  4  and  5 . 75 

Figure  30.  WT-based  denoising  for  constant  amplitude  chirp  type,  4  level,  method  1, 

Daubechies  wavelet  4,  SNR=-6  dB . 76 

Figure  31.  WT-based  denoising  for  constant  amplitude  chirp  type,  4  level,  method  1, 

Daubechies  wavelet  4,  SNR=0  dB . 77 

Figure  32.  WT-based  denoising  for  constant  amplitude  chirp  type,  4  level,  method  1, 

Daubechies  wavelet  4,  SNR=6  dB . 78 

Figure  33.  WT-based  denoising  for  increasing  amplitude  chirp  type,  4  level,  method 

1,  with  Daubechies  wavelet  orders  3,  4  and  5 . 79 

Figure  34.  WT-based  denoising  for  increasing  amplitude  chirp  type,  4  level,  method 

1,  Daubechies  wavelet  5,  SNR=-6  dB . 80 

Figure  35.  WT-based  denoising  for  increasing  amplitude  chirp  type,  4  level,  method 

1,  Daubechies  wavelet  5,  SNR=0  dB . 81 

Figure  36.  WT-based  denoising  for  increasing  amplitude  chirp  type,  4  level,  method 

1,  Daubechies  wavelet  5,  SNR=6  dB . 82 

Figure  37.  Wavelet  shrinkage  algorithm  and  proposed  WT-based  denoising  scheme 

performances  for  the  sinusoidal  signal  type,  using  3 -level  decomposition 

and  Daubechies  wavelet  5 . 85 

Figure  38.  Wavelet  shrinkage  algorithm  and  proposed  WT-based  denoising  scheme 
performances  for  the  constant  amplitude  chirp  type,  using  4-level 

decomposition  and  Daubechies  wavelet  4 . 86 

Figure  39.  Wavelet  shrinkage  algorithm  and  proposed  WT-based  denoising  scheme 

performances  for  the  increasing  amplitude  chirp  type,  using  4-level 

decomposition  and  Daubechies  wavelet  5 . 87 


xii 


ACKNOWLEDGMENTS 


I  wish  to  express  my  sincerest  appreciation  to  my  advisors  Professor  Monique 
Fargues  and  Professor  Ralph  Hippenstiel  for  their  support  and  supervision. 

I  would  like  to  thank  Professor  Roberto  Cristi  for  his  contributions. 

Most  importantly,  I  would  like  to  thank  my  wife  for  her  love,  care  and  support.  I 
would  not  have  been  able  to  do  this  without  her. 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


xiv 


EXECUTIVE  SUMMARY 


Hie  main  goal  of  this  thesis  was  to  develop  a  denoising  scheme  to  recover  signals 
that  are  distorted  by  additive  white  Gaussian  noise,  regardless  of  the  signal’s  sperctral 
content.  Two  approaches  were  investigated,  a  short-time  Fourier  transform  and  a  wavelet 
transform-based  approach.  Both  transforms  are  used  to  decompose  a  signal  into  several 
frequency  bands  and  to  minimize  the  noise  in  each  band.  Each  frequency  band  content  is 
tested  for  Gaussianity,  which  is  accomplished  by  investigating  the  signal’s  kurtosis.  It  is 
well  known  that  the  signal  kurtosis  is  equal  to  three  only  when  the  signal  is  Gaussian  with 
zero  mean.  Therefore  this  parameter  can  help  to  provide  a  reliable  determination  whether 
or  not  a  data  sequence  is  Gaussian.  However,  the  kurtosis  estimate  of  a  short  data 
segment  is  not  always  reliable.  The  Bootstrap  method  is  used  in  the  kurtosis  estimation  to 
overcome  this  difficulty.  The  Bootstrap  is  a  statistical  scheme,  which  improves  the 
reliability  of  a  parameter  estimate  in  situations  where  conventional  techniques  are  not 
valid  because  of  short  data  length  issues.  The  Bootstrap  uses  data  sequences  that  have  the 
same  length  as  the  original  signal,  but  are  obtained  by  randomly  resampling  the  data 
using  replacement.  These  resampling  is  done  many  times  and  each  data  set  is  treated  as 
repeated  experiment.  Then,  the  parameter  of  interest  is  estimated  for  each  of  these 
resampled  sequences,  to  obtain  a  statistic  for  the  parameter  of  interest.  This  resulting 
statistic  can  then  be  used  to  perform  hypothesis  tests  on  the  parameter  value. 

If  the  band  is  estimated  to  be  Gaussian  (i.e.,  the  noise  is  dominant),  then 
thresholding  is  applied  to  the  band  to  minimize  the  noise  effects.  Two  thresholding 
techniques  are  considered,  hard  and  soft  threshold.  Hard  thresholding  coefficients  are  0 
for  Gaussian  data  and  1  for  non- Gaussian  data.  Whereas  soft  thresholding  coefficients  are 
obtained  by  considering  the  data’s  closeness  to  Gaussianity.  The  closer  to  being  Gaussian 
the  band  specific  data  is,  the  smaller  the  soft  thresholding  coefficients  are  and  vice  versa. 
The  denoised  signal  is  obtained  by  applying  the  appropriate  inverse  transformation. 

Three  different  test  signal  types  are  selected  to  investigate  the  performances:  a 
sinusoid,  a  chirp  with  constant  amplitude,  and  a  chirp  with  an  RC  time  constant -like 


xv 


amplitude  increase.  The  mean  square  error  (MSE)  and  a  distance  measure  defined 
between  original  and  recovered  signals  are  selected  to  compare  performances.  Results 
show  FFT-based  denoising  schemes  perform  better  than  WT-based  denoising  schemes  on 
the  stationary  sinusoid  signal  type,  whereas  WT-based  schemes  outperform  FFT-based 
schemes  on  chirp  type  signals.  Finally,  results  show  that  the  soft  thresholding  scheme 
always  performs  at  least  as  well  as  or  better  than  the  hard  thresholding  one. 


I.  INTRODUCTION 


In  many  data  transmission  and  storage  systems,  noise  gets  introduced  into  data, 
which  reduces  the  signal  quality.  As  a  result,  various  filtering  techniques  have  been 
designed  to  suppress  noise  contributions  in  order  to  improve  the  overall  signal  quality. 

Fourier  and  Wavelet  transforms  decompose  a  noisy  signal  into  several  frequency 
bands.  Traditional  filter  design  methods  have  requirements  on  the  frequency,  magnitude 
and  phase  of  the  signal  such  as  passband  ripple,  stopband  attenuation,  transition  width, 
and  phase  constraints.  The  assumption  behind  these  design  criteria  is  that  the  signal  is 
restricted  to  be  in  a  certain  frequency  band  and  that  the  frequencies  outside  this  specific 
band  are  treated  as  distortion.  Note  that  this  paradigm  breaks  down  when  signal  and 
distortion  terms  overlap  in  frequency. 

Denoising  attempts  to  remove  the  noise  and  to  recover  the  original  signal 
regardless  of  the  signal’s  frequency  content.  The  basic  idea  is  to  look  at  each  frequency 
band  of  interest  and  to  minimize  its  noise  effect  by  retaining  the  dominant  component. 

The  band  is  left  untouched  when  the  signal  is  dominant  so  as  not  to  lose  the  signal 
component,  while  thresholding  is  applied  when  the  noise  is  dominant. 

This  thesis  discusses  a  denoising  scheme  that  implements  frequency  band  specific 
thresholding  schemes  using  the  Bootstrap  method  and  the  kurtosis.  Chapter  2  presents  the 
processing  techniques  and  the  Bootstrap  method.  Chapter  3  discusses  the  proposed 
denoising  scheme.  Simulation  results  are  presented  in  Chapter  4.  Finally,  Chapter  5 
presents  conclusions. 
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H.  BACKGROUND 


A.  PROCESSING  TECHNIQUES 

Signal  processing  allows  users  to  extract  relevant  information  from  a  given  signal. 

When  the  raw  data  does  not  allow  extracting  the  desired  information,  transformations  to 
another  domain  may  be  performed  to  do  so.  The  most  common  transformation  types  are 
the  Fourier  and  Wavelet  transforms,  which  are  discussed  in  the  following  two  sections. 

1.  Fourier  Analysis 

Fourier  analysis  allows  the  representation  of  a  given  signal  as  a  linear 
combination  of  complex  sinusoids  with  different  frequencies.  For  periodic  signals  this  is 
called  the  Fourier  Series.  This  representation  becomes  extended  to  the  Fourier 
Transform  when  the  signal  is  aperiodic.  Both  representations  are  discussed  in  the 
following  subsections. 

a  Fourier  Series 

A  periodic  signal  x(t)  with  period  T0  may  be  represented  as  an  infinite 
linear  combination  of  complex  exponentials: 

jk—t 

x(t)=Y,ake  T°  .  (2.1) 

k  =— oo 

271 

jk — t 

where  /0  =  I  /  T{)  is  One,  fundamental  frequency,  e  is  called  the  k'h  harmonic,  and  ak  is 
the  k,h  weight.  This  representation  is  called  the  Fourier  series  representation  [1],  where 
the  set  of  coefficients  ak  are  called  the  Fourier  series  coefficients  which  are  obtained  by 

1  r  -7*— t 

ak  =  —  J  x(t)e  T°  dt,  k  =  — (2.2) 

T0  T0 

The  coefficients  af  s  are  a  measure  of  the  strength  of  the  signal’s  components  at  the  k,b 
harmonic  of  the  fundamental  frequency. 
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b.  Fourier  Transform 

For  aperiodic  signals  the  Fourier  series  representation  can  be  extended. 

Assuming  that  xp  (t)  is  a  periodic  signal  with  a  period  T0  ,  and  x(t)  represents  one  period 

of  xp  (t)  ,  then  as  T0  increases  xp  (t)  is  identical  to  x(t)  over  a  longer  interval.  Therefore, 
in  the  limiting  case  [2] 

x(t)  =  lim  x  (t) .  (2.3) 

Tp^oo 

Replacing  the  limits  of  the  integral  in  equation  (2.2)  by  and  °° ,  and  multiplying  both 

sides  by  Tp  ,  leads  to  the  Fourier  Transform  of  x( t)  given  by 

X(f)=  ]  x(t)e~j2%fidt,  (2.4) 


where  f  -  k/Tp  and  X{f)-Tpak.  The  Inverse  Fourier  Transform  is  used  to  recover  the 
original  signal  x(t)  from  X(f) ,  and  is  defined  as 

x(t)  =  j  X(f)ei2%ftdf  .  (2.5) 


signals  as 


Hie  Discrete  Time  Fourier  Transform  (DTFT)  is  defined  for  discrete- time 


X(a>)  -^x(n)  e  , 

n—oo 


(2.6) 


and  the  corresponding  inverse  is  given  by 


x(n)  = —  [  X (CO)  ejm"cl(j)  . 
2n  J 
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(2.7) 


It  should  be  noted  that  the  DTFT  is  periodic  with  period  2k  . 

The  Discrete  Fourier  Transform  (DFT)  for  a  finite  time  discrete-time 
signal  x( n )  with  n=0,...,  N  - 1  is  defined  as 

X(k)  =  Y^xin)  e-j2nkn,N  .  (2.8) 

n—0 


4 


Its  inverse  is  given  by 


x(n)  =  ^YiX(k)ei*lhllN.  (2.9) 

The  Fast  Fourier  Transform  is  a  computationally  efficient  implementation  of  the 
Discrete  Fourier  Transform  when  the  signal  length  N  is  a  power  of  two.  The  Discrete 

'j 

Fourier  Transform  requires  N  multiphcations  whereas  the  Fast  Fourier  Transform 
requires  (iV72)log2  N  multiplications. 

c.  Short-Time  Fourier  Transform 

Hie  basic  Fourier  Transform  allows  the  passage  from  the  time  domain  to 
the  frequency  domain.  However,  one  may  need  to  preserve  both  the  time  and  the 
frequency  information  contained  in  non- stationary  signals.  Unfortunately,  the  basic 
Fourier  Transform  does  not  provide  such  dual  information. 

Time  localization  can  be  introduced  by  windowing  the  signal  before  using 
the  Fourier  Transform,  where  the  window  size  is  selected  short  enough  to  assume  the 
signal  inside  the  window  is  stationary.  Then,  the  Fourier  transform  may  be  implemented 
on  each  windowed  signal  portion.  The  resulting  transformation  is  called  the  Short-Time 
Fourier  Transform  and  expressed  as 

S(x ,  / )  =  J x(t)w(t  -x )  e~i2nft dt ,  (2.10) 

t 

where  w(t  — x )  denotes  the  sliding  window  centered  around  x  .Note  that  S(% ,  f) 
provides  a  two-dimensional  representation  of  the  signal  frequency  information  at  various 
times  x  and  is  a  “local”  spectrum  of  the  signal  x(t)  around  the  analysis  point  x  [3]. 

Many  different  window  types  may  be  selected,  depending  on  the  characteristics  of  the 
Short-Time  Fourier  Transform  desired.  In  addition,  the  window  length  determines  the 
resolution  in  time  and  in  frequency;  good  time  resolution  requires  a  short  window, 
whereas  good  frequency  resolution  requires  a  long  window.  The  joint  time- frequency 
resolution  of  the  Short-Time  Fourier  Transform  is  limited  by  the  uncertainty  principle.  A 
short  time  window  results  in  a  loss  of  frequency  resolution,  and  vice  versa. 
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2.  Wavelet  Analysis 

A  wavelet  is  defined  as  a  small  wave  which  has  its  energy  concentrated  in  time 
and  frequency  [4].  Such  a  characteristic  is  useful  for  the  analysis  of  transient,  non¬ 
stationary  or  time -varying  phenomena.  By  comparison,  sinusoidal  functions  used  in 
Fourier  Analysis  have  a  constant  amplitude. 

The  Wavelet  Transform  (WT)  provides  an  alternative  to  the  Short-Time  Fourier 
Transform  (STFT),  as  it  uses  short  windows  at  high  frequencies  and  long  windows  at  low 
frequencies  [5],  while  the  STFT  uses  windows  of  constant  size,  as  illustrated  in  Figure  1. 
Note  that  the  Wavelet  Transform  still  satisfies  the  uncertainty  principle,  however,  the 
time  resolution  becomes  arbitrarily  good  at  high  frequencies,  while  the  frequency 
resolution  becomes  arbitrarily  good  at  low  frequencies.  The  following  subsections 
discuss  different  types  of  Wavelet  Transform  that  are  used  for  continuous  or  discrete  time 
signals. 


frequency  frequency 


(a)  STFT  (b)  WT 


Figure  1.  (a)  Time -Frequency  resolution  for  the  STFT,  (b)  Time -Frequency  resolution  for  the  WT. 
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a.  The  Continuous  Time  Wavelet  Transform 

The  Continuous  Time  Wavelet  Transform  may  be  applied  to  continuous 
signals  and  two  dimensions  in  the  transfer  domain  are  continuous.  It  is  similar  to  the 
Fourier  Transform  in  that,  it  is  obtained  by  projecting  the  signal  onto  a  basis  function. 
However,  the  Continuous  Time  Wavelet  Transform  projects  the  signal  onto  scaled  and 
shifted  versions  of  the  wavelet  function  while  the  Fourier  Transform  uses  complex 
exponentials  as  basis  functions.  The  Continuous  Time  Wavelet  Transform  is  defined  as 

C0,t)  =  4=  ?  x(t)\\r(— -)dt, 
sjs  s 

where  \\t  (t)  is  the  wavelet  function,  x  denotes  the  shift  in  time  and  s  is  the  scale  factor 
that  denotes  compression  or  expansion  in  time.  The  inverse  transform  for  finite  K  is 
obtained  by 

x(t)  =  K  f  f  -\c(s,T  )\|/  (- — —)dsdx  , 

J  J  s  s 

where  the  parameter  K  is  given  by 


and  T  to )  is  the  Fourier  Transform  of  the  wavelet  function  \|/  (t) .  [6] 
b.  The  Discrete  Time  Wavelet  Transform 

The  Discrete  Time  Wavelet  Transform  is  the  discrete  time  version  of  the 
Continuous  Wavelet  Transform.  It  is  used  for  discrete  time  signals  and  the  dimensions  of 
the  transform  domain  are  discrete  as  well.  The  Discrete  Wavelet  Transform  is  obtained  as 

C(a,b)  =  V-^x(n)\|/ (— -) , 

„  si  a  a 

where  a ,  b  and  n  are  the  discrete  parameter  versions  of  s,  X  and  1  given  in  Equation 
(2.1 1),  respectively.  The  scaling  factor  a  has  another  restriction  as  a-aJQ  where 


(2.11) 


(2.12) 


(2.13) 


(2.14) 
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7=0,1,. .  ,,log2 (signal  length).  The  common  choice  for  a0  is  2,  as  it  allows  for  fast 
algorithms. 

c.  Mallat’s  Algorithm 

The  Discrete  Wavelet  Transform,  can  be  implemented  by  using  Mallat’s 
algorithm  when  a0  =  2  [4],  and  is  illustrated  in  Figure  2  for  a  three- level  decomposition. 
The  signal  is  passed  through  a  high-pass  and  a  low-pass  filter,  both  of  which  have  a 
bandwidth  of  half  the  signal  spectrum.  Then,  following  the  Nyquist’s  rule,  the  outputs  of 
the  filters  are  subsampled  by  two.  The  subsampled  high- pass  filter  output  is  called  the 
Detail  Sequence  and  the  subsampled  low- pass  filter  output  is  called  the  Approximation 
Sequence.  This  procedure  may  be  recursively  applied  to  the  approximation  sequence 
obtained  at  previous  levels.  Note  that  a  signal  of  length  2  '  can  be  decomposed  only  j 
times,  because  the  approximation  sequence  has  only  one  sample  left  after  j  levels. 


Figure  2.  Mallat’s  algorithm. 


B.  HIGHER  ORDER  STATISTICS 
1.  Moments  and  Cumulants 

The  first  four  moments  for  a  real  valued  and  stationary  signal  x(n)  are  given  by 


m,  =  E{x(n )} , 

(2.15) 

! )  =  E{x  (n  (n  +T1 )} , 

(2.16) 
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m3(x1  ,x2)  =  E{x(ri)x(n  +x1)x(n+x2)] , 


(2.17) 


m  4  (x  j  ,x  2  ,x  3 )  =  E{x(ri)x(n  +x,)x(n  +x2)x(n  +X3)] .  (2.18) 

The  first  two  moments  are  equal  to  the  mean  and  correlation  functions  respectively.  The 
first  four  cumulants  are  given  by 

=  m, .  (2.19) 

c2 (x, )  =  m2( t ,)  - ml (Xj)2,  (2.20) 

c3(x,  ,x2)  =  m3(x1  ,x2)-ml[m1  (x,)  +  m2(x2)  +  m2(x2  -x3)]  +  2(m,)3,  (2.21) 

c4(Xj  ,x2,x3)  m4(x1,x2,x3)  m2(x,)m,(x3  x2)  /??2 (x 2)/?r  (x3  Xj) 

-  m2 (x  3)m,  (x2  -x  J  -  m,  [m,  (x2  -x ,  ,x3  -xt ) 

+  m3 (x px2 )  +  m3  (x,  ,x3 )  +  m3 (x2  ,x3 )]  (2.22) 

+  (m,)2[m2(x1)  +  m2(x2)  +  m2(x3)  +  m2(x3  -xt) 

+  m2  (x3  — X2 )  +m2(x2  -X, )]  -  b^!)4. 

Note  that  the  second  and  third  order  cumulants  are  identical  to  the  second  and  third  order 
moments  respectively  when  x(n)  is  a  zero-mean  process. 

Cumulants  have  properties  that  make  them  more  desirable  than  moments.  Some 
of  these  properties  are 

a.  Each  cumulant  is  independent  of  all  lower  order  cumulants. 

b.  All  cumulants  of  order  greater  than  two  are  equal  to  zero  for  Gaussian 
processes.  Hence,  any  Gaussian  process  is  completely  characterized  by  its  first  two 
cumulants.  Therefore  higher- order  cumulants  can  be  used  to  estimate  the  degree  of  non- 
Gaussianity  of  a  process. 

c.  Cumulants  of  the  sum  of  two  independent  statistical  processes  are  equal  to 
the  sum  of  their  respective  cumulants.  [6] 

2.  Variance,  Skewness  and  Kurtosis  Measures 

Setting  Xj,x2,x3  equal  to  0  in  (2.20),  (2.21),  (2.22),  and  assuming  that  mj= 0,  leads 

to  the  variance  yl ,  skewness  y2 ,  and  kurtosis  y3  measures: 
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y2  =  E{x(n)2}  =c2(0), 


(2.23) 


_E{x(nf  j 

Y3  _  3  ’ 

o 

E{x(n)4  j 

y4=  4  . 


where  a  is  the  standard  deviation  around  the  mean  of  the  signal.  [7] 


(2.24) 

(2.25) 


C.  THE  BOOTSTRAP 

The  Bootstrap  is  a  powerful  technique  for  assessing  the  accuracy  of  a  parameter 
estimator  in  situations  where  conventional  techniques  are  not  applicable  [8].  In  many 
applications,  one  needs  to  estimate  one  or  more  parameters  of  a  random  process,  and/or 
calculate  some  statistical  parameters  such  as  the  mean  or  variance.  Most  of  the  estimation 
techniques  used  for  this  purpose  assume  that  the  set  of  samples  used  in  the  estimation  is 
large  enough  to  reach  asymptotic  results.  However,  in  practice  this  assumption  usually 
does  not  hold  as  the  sample  set  may  not  be  large  enough  or  the  samples  may  be  non¬ 
stationary.  The  Bootstrap  scheme  randomly  reassigns  the  observations,  recomputes  the 
estimates  many  times,  and  treats  these  reassignments  as  repeated  experiments.  In  this 
section,  first  the  basic  Bootstrap  principle  is  stated.  Next,  the  usage  of  Bootstrap  in 
estimating  the  confidence  interval  for  a  parameter  is  discussed,  and  finally  this  discussion 
is  extended  to  hypothesis  testing. 

1.  Basic  Principle 

Let  x-{x  1,  A' 2,  •  •  • ,  xnJ  be  a  collection  of  n  independent  and  identically  distributed 
random  variables  drawn  from  an  unknown  distribution,  D.  Let  p  denote  an  unknown 

statistical  parameter  of  D  such  as  the  mean  or  the  variance,  and  p  denote  an  estimator  of 
p,  calculated  from  a.  If  the  estimate  p  is  to  be  used  in  place  of  the  real  parameter  p,  it 

may  be  important  to  know  the  sampling  distribution  of  p  .  The  distribution  of  p  may  be 
estimated  with  the  Bootstrap  method,  and  is  obtained  by  resampling  many  times  from  a 

distribution  D ,  chosen  to  be  close  to  D.  such  that  I)  approaches  D  as  ft  — >  .  Note  that 
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the  choice  of  I)  is  not  unique.  If  the  type  of  D  is  known  but  its  statistical  parameters  are 

not  known,  then  I)  is  chosen  as  a  distribution  of  the  same  type  as  D ,  with  the  statistical 
parameters  obtained  from  x.  For  example,  if  we  know  that  the  data  is  Gaussian  but  do  not 
know  its  mean  and  variance,  we  perform  the  resamplings  assuming  the  data  to  have  the 
same  mean  and  variance  as  x  has.  This  approach  is  called  the  parametric  Bootstrap.  If 
nothing  is  known  about  D,  the  resamplings  are  drawn  from  x  with  replacement,  so  that 
each  value  in  a  resample  set  has  probability  equal  to  1/n.  This  approach  is  called  the 
nonparametric  Bootstrap. 

2.  Parameter  Confidence  Interval 

The  Bootstrap  principle  may  also  be  applied  to  obtain  a  (l-a)100%  confidence 
interval  for  the  parameter  p.  First,  the  data  set  is  resampled  many  times  such  that  each 
resampled  set  is  of  the  same  size  as  the  data.  Next,  an  estimator  pk  is  obtained  from  each 
resampled  set,  where  k  -  1 ....,  N  and  N  is  the  number  of  repetitions.  Then,  the  estimates 

pk  are  sorted  in  increasing  order.  Finally  the  indices  of  the  lower  limit  pL  and  the  upper 
limit  p,,  of  the  estimates  in 

P(PL  ^  P^Pu)  =  l~a 


are  obtained  by  using 


L  = 


Na 


and 


U=N-  L+  1, 


where  |_Aj  denotes  the  integer  part  of  the  value  A. [8] 

3.  Hypothesis  Testing 
a.  Description 

Suppose  one  needs  to  perform  a  hypothesis  test,  such  as  H  :  p  <  p0 
against  the  hypothesis  K  :  p>  p0,  where  p0  is  given.  A  new  statistic,  defined  as 


(2.26) 


(2.27) 


(2.28) 


11 


(2.29) 


j-P-J  o 

o 

may  be  selected,  where  a  is  the  estimator  of  the  standard  deviation  a  of  p  .  The 
estimator  for  the  standard  deviation  o  of  p  will  be  defined  later. 

To  perform  the  hypothesis  testing,  one  first  draws  random  sequences 
X,  ,  x*2 ,...,  xN,  of  the  same  size  as  x,  with  replacement  from  x.  Note  that  4  does  not  denote 
complex  conjugation,  but  means  that  this  is  a  resampled  set,  or  a  parameter  obtained  from 

xv*  ^ 

a  resampled  set.  Then  the  statistic  T  is  estimated  from  each  sequence  x  as 

T  (2.30) 

o 

where  p  and  a  are  estimated  parameters  obtained  from  the  resample  x* ,  instead  of  x, 
and  the  constant  pQ  is  replaced  with  p  .  Note  that  o  is  included  as  a  scale  factor  in  the 

calculation  of  T  .  Dividing  by  o  is  called  Bootstrap  pivoting  and  it  is  done  to  ensure  T 
is  asymptotically  pivotal  when  n  — >  °° ,  i.e.,  the  asymptotic  distribution  of  T  does  not 
depend  on  any  unknown  parameters.  Replacing  p0  with  p  and  using  Bootstrap  pivoting 

yy  *  XV  *  XV  XV  * 

is  important  because  the  Bootstrap  distribution  of  T  =(p  -  p)/o  is  a  better 
approximation  to  the  distribution  of  T  =  (p-  p0)/ a  under  H ,  than  the  Bootstrap 

XV*  XV  XV 

distribution  of  S*  -  p  -p  is  to  the  distribution  of  S  =  p  —p0  under  H  [8,  9].  Next,  the 

XV*  XV  *  XV  * 

set  of  test  statistics  Ti  N  are  sorted  by  increasing  order,  and  the  hypothesis  H  is 

XV  XV  * 

rejected  if  T  where  M  is  chosen  according  to  N  and  the  level  of  significance  a  as 

[8,  pp.  62] 

M  =  (N  +l)(l-a) .  (2.31) 
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Note  that  the  test  statistic  T  is  given  by 


P  ~  P  o 
o 


(2.32) 


when  the  hypotheses  to  be  tested  are  H  :  p  =  p0  against  K  :  p  ^  p0. 

b.  Estimation  of  the  Standard  Deviation  for  p 

The  parameter  o  can  be  estimated  by  using  the  Bootstrap.  Towards  that 
end,  resamples  of  the  same  length  as  the  data  set  x,  called  x*  are  drawn  from  x  randomly 
with  replacement,  to  obtain  a  total  of  B  resamples.  After  resampling,  the  Bootstrap 

XV  *  XV 

estimates  p  are  calculated  in  the  same  manner  as  p  was,  but  with  resamples  x*  instead 
of  x.  As  a  result,  the  standard  deviation  o  of  p  is  estimated  by 

n  l:  '  (2'33) 

In  this  chapter  we  presented  the  processing  techniques  that  we  used  in  the 
proposed  denoising  scheme.  In  the  next  chapter  we  will  discuss  the  proposed  denoising 
scheme. 
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HI.  DENOISING  USING  THE  KURTOSIS  AND  THE  BOOTSTRAP 


A.  INTRODUCTION 

Fourier  and  Wavelet  transforms  decompose  a  signal  into  several  frequency  bands. 
As  a  result,  white  Gaussian  noise  components  affect  all  frequency  bands.  The  proposed 
technique  examines  each  frequency  band  and  tends  to  minimize  the  white  noise 
contribution.  First,  we  obtain  the  signal  time-frequency  representation  (or  time- scale 
representation  in  the  Wavelet  transform  case).  Next,  each  frequency  band  is  tested  for 
Gaussianity,  and  thresholding  is  performed  on  the  spectral  location  found  to  be  Gaussian. 
Last,  the  signal  is  transformed  to  the  time  domain  using  the  appropriate  inverse 
transform.  The  following  sections  discuss  two  different  implementations  using  Fourier 
and  Wavelet  transformations. 

B.  FAST  FOURIER  TRANSFORM  BASED  DENOISING 

The  FFT-based  denoising  scheme  is  performed  in  four  steps,  as  illustrated  in 
Figure  3.  For  computational  convenience,  the  signal  length  is  assumed  to  be  512  points. 

A  longer  data  set  could  have  been  used,  resulting  in  more  frequency  bands,  and  more 
computations.  A  data  set  of  minimum  length  of  512  allows  separation  in  time  and 
frequency.  A  shorter  segment  will  not  provide  a  reasonable  number  of  time -frequency 
cells. 


Figure  3.  FFT-based  denoising  scheme. 
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1.  Short-Time  Fourier  Transform  Step 

In  the  first  step,  the  data  is  divided  into  32  data- point  segments  and  weighted  by  a 
triangular  window.  A  four- to- one  overlap  is  used  to  obtain  a  reasonable  data  size.  Then, 
each  segment  is  transformed  to  the  frequency  domain  with  the  Fast  Fourier  Transform. 

This  results  in  17  frequency  bands.  A  larger  data  set  would  have  allowed  for  a  larger 
segment  size.  The  STFT  information  is  contained  in  a  matrix  of  dimensions  61  by  32 
elements  for  a  data  length  of  5 12  points,  where  the  first  number  corresponds  to  the  time 
dimension  and  the  second  one  corresponds  to  the  frequency  dimension. 

2.  Gaussianity  Test  Step 

The  second  step  tests  the  Gaussianity  of  the  transformed  values  in  each  frequency 
band,  where  real  and  imaginary  parts  of  the  data  are  tested  separately  with  the  kurtosis. 
Recall  that  the  normalized  kurtosis  value  for  a  Gaussian  data  sequence  is  equal  to  3. 

Hence,  looking  at  the  kurtosis  value  may  give  an  idea  of  a  sequence’s  Gaussianity, 
provided  the  data  length  is  sufficient  for  the  estimation  to  be  meaningful.  However,  the 
sequences  in  this  particular  scheme  are  of  length  equal  to  61,  which  may  not  be  sufficient 
for  meaningful  estimations.  Therefore,  the  Gaussianity  test  is  implemented  as  a 
hypothesis  test,  with  the  hypothesis  H  \kurtosisof  thedata  =  3  ,  and  the  hypothesis 

K  \kurtosisof  thedata  ^  3  .  The  sequence  is  found  to  be  Gaussian  when  the  hypothesis  H 

is  true,  meaning  the  kurtosis  value  for  the  data  sequence  is  3  within  the  specified 
confidence  interval.  However,  the  sequence  is  found  to  be  non- Gaussian  when  hypothesis 
K  is  true,  meaning  the  kurtosis  is  not  equal  to  3  within  the  selected  confidence  interval. 

This  hypothesis  testing  is  performed  using  the  Bootstrap  method  discussed  earlier  in 
Chapter  2.  An  empirical  confidence  interval  a  =  0.05  was  selected  here. 

3.  Thresholding  Step 

Signal  components  are  found  in  a  particular  sequence  when  some  frequency 
bands  are  estimated  to  be  non- Gaussian  for  both  real  and  imaginary  parts.  Thresholding  is 
applied  to  each  band  estimated  to  be  Gaussian,  to  minimize  noise  effects  in  these  bands. 
Two  thresholding  schemes  are  considered  below. 
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a. 


Hard  Thresholding 


Hai'd  thresholding  zeroes  out  all  values  in  the  frequency  band  that  is  found 

to  be  Gaussian.  The  hard  thresholding  coefficient  is 

f  0,  if  thebandisGaussian 
ch=\ 

1,  if  thebandisnotGaussian. 


(3-1) 


b.  Soft  Thresholding 


Soft  Thresholding  is  obtained  by  multiplying  values  in  the  specific 
frequency  band  that  is  found  to  be  Gaussian,  by  a  coefficient  between  0  and  1.  Using  a 
coefficient  of  0  is  the  same  as  hard  thresholding,  whereas  using  a  coefficient  of  1  is  the 
same  as  leaving  the  frequency  band  undisturbed.  The  soft  thresholding  coefficient  is 
calculated  using 


(3.2) 


where  y4  is  the  bootstrapped  kurtosis  of  the  particular  frequency  band,  and  |  |  denotes 

the  absolute  value.  The  bootstrapped  kurtosis  value  y4  is  limited  not  to  exceed  4.5,  which 

will  be  explained  later.  The  bootstrapped  kurtosis  value  y4  is  obtained  by  using  the 
Bootstrap  principle,  which  calls  for  resamplings  from  the  data  set  many  times  with 
replacement,  to  obtain  N  resampled  sets  of  the  same  length  as  the  original  data  set.  Next, 

the  kurtosis  value  for  each  resample  is  found.  Finally,  the  Bootstrapped  kurtosis  y4  is 
defined  as  the  estimated  mean  obtained  from  N  kurtosis  values.  It  should  be  noted  that  the 
thresholding  coefficient  c  is  a  function  of  the  frequency  band’s  degree  of  Gaussianity. 
Equation  (3.1)  shows  that  the  coefficient  c  gets  closer  to  0  as  the  bootstrapped  kurtosis 

value  y4  for  a  specific  frequency  band  gets  closer  to  the  theoretical  value  3,  and  vice 
versa.  Therefore,  the  closer  a  frequency  band  gets  to  being  Gaussian,  the  smaller 
contribution  it  has  after  soft  thresholding. 

Recall  that  frequency  bands  to  be  thresholded  are  those  that  passed  the 
Gaussianity  test.  Therefore,  one  would  expect  their  corresponding  bootstrapped  kurtosis 
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estimates  to  be  close  to  the  theoretical  value  3.  However  in  some  cases,  the  estimated 
bootstrapped  kurtosis  value  may  be  far  off  from  3.  For  example,  a  few  values  in  the  band 
may  be  higher  than  the  others  especially  when  a  frequency  bin  contains  a  short  data 
segment.  In  this  case,  the  band  may  pass  the  Gaussianity  test  as  most  of  the  band  is 
Gaussian,  but  may  still  have  a  bootstrapped  kurtosis  value  high  enough  to  obtain  a 
thresholding  coefficient  c  greater  than  1.  Note  that  a  thresholding  coefficient  greater  than 
1  would  amplify  noise  contributions  in  the  frequency  band,  instead  of  suppressing  them. 

Therefore,  the  bootstrapped  kurtosis  value  y4  is  limited  to  4.5  in  (3.1)  to  avoid  such 
potential  noise  amplification.  This  value  is  chosen  empirically.  Thus,  if  the  bootstrapped 
kurtosis  value  of  a  specific  frequency  band  is  greater  than  4.5,  the  data  point  with  the 
largest  absolute  value  in  that  band  is  stored  away  and  replaced  with  a  zero.  Then,  the 
bootstrapped  kurtosis  for  that  band  is  estimated  again.  This  procedure  is  repeated  until  a 
kurtosis  value  less  than  or  equal  to  4.5  is  obtained.  However,  the  number  of  repetitions 
should  be  limited  to  ensure  that  most  of  the  data  is  not  zeroed  out,  to  prevent 
bootstrapped  kurtosis  estimation  problems.  An  empirical  limit  of  one  third  of  the  number 
of  time  frames  defined  in  the  STFT.  Note  that  if  the  estimated  kurtosis  value  is  still 
greater  than  4.5  after  that  many  iterations,  then  it  is  set  to  4.5,  which  provides  a 
thresholding  coefficient  equal  to  1.  This  causes  no  change  on  the  frequency  band  under 
consideration. 

After  multiplying  the  frequency  band  with  the  thresholding  coefficient,  the 
values  that  were  removed  to  limit  the  kurtosis  estimation  are  reinserted  into  their  original 
locations.  This  allows  for  the  non- Gaussian  values,  which  may  correspond  to  the  signal 
components,  to  not  be  affected  by  the  thresholding  step. 

4.  Inverse  Fourier  Transform 

Real  and  imaginary  parts  of  the  frequency  bands  are  combined  to  form  the  time- 
frequency  representation  matrix  after  the  thresholding  step,  and  the  recovered  signal  is 
obtained  using  the  inverse  Fourier  transform. 
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C.  WAVELET  TRANSFORM  BASED  DENOISING 

The  Wavelet  transform- based  denoising  scheme  is  similar  to  the  FFT-based 
scheme,  as  illustrated  in  Figure  4.  Note  that  only  real  valued  transforms  are  obtained  in 
this  case  due  to  the  real  wavelet  transform  kernel  form  selected  in  our  work  (Daubechies 
wavelets  of  order  3,  4  or  5).  The  choice  of  Daubechies  wavelets  was  motivated  by  the 
earlier  works  in  GSM  signal  denoising  reported  by  Aktas  and  Mantis  [10,  11].  This 
denoising  scheme  is  performed  in  four  steps. 


Figure  4.  WT-based  denoising  scheme. 


1.  The  Wavelet  Transform 

First,  the  Wavelet  transform  is  applied  to  the  noisy  signal.  Three-,  four-,  or  five- 
level  decompositions  are  considered  in  this  work,  because  simulations  showed  that 
higher- level  decompositions  do  not  improve  performance,  given  the  data  length 
considered  (512). 

2.  Gaussianity  Test 

The  approximation  and  detail  coefficients  of  a  Gaussian  data  set  remain  Gaussian 
[12].  Thus,  detail  and  approximation  coefficients  are  tested  for  Gaussianity  and  a  decision 
is  made  for  each. 
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3.  Thresholding 

The  same  thresholding  procedure  as  that  was  considered  earlier  for  the  FFT-based 
denoising  scheme  is  applied  to  the  detail  coefficients.  However,  thresholding  the 
approximation  coefficients  may  result  in  removing  a  significant  portion  of  the  signal.  One 
can  use  one  of  the  following  four  schemes  to  threshold  the  approximation  coefficients  to 
minimize  potential  distortion:  Apply  the  thresholding  scheme  as  that  used  on  the  detail 
coefficients,  leave  coefficients  unperturbed,  filter  the  approximation  coefficients  using  a 
median  filter  of  order  three,  or  use  a  predictor  of  order  two.  The  last  two  options  were 
investigated  as  ways  to  smooth  the  approximation  coefficients  but  they  did  not  perform 
as  good  as  the  first  two.  Results  for  these  operations  are  discussed  in  the  next  chapter. 

4.  Inverse  Wavelet  Transform 

Finally,  updated  approximation  and  detail  coefficients  are  inverse  Wavelet 
transformed  to  obtain  the  recovered  signal. 

In  this  chapter  the  proposed  denoising  scheme  was  discussed  in  detail.  This 
scheme  was  tested  using  three  test  signals.  The  test  signal  descriptions  and  the  simulation 
results  are  presented  in  the  next  chapter. 
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IV.  SIGNAL  DESCRIPTION  AND  SIMULATION  RESULTS 


The  denoising  scheme  presented  in  the  previous  chapter  was  tested  using  three 
test  signals.  The  codes  used  for  the  simulations  are  presented  in  Appendix  A.  In  this 
chapter  the  test  signals  are  described  and  the  simulation  results  are  summarized. 

A.  SIGNAL  DESCRIPTION 

Three  different  types  of  test  signals  are  used:  a  sinusoid,  a  chirp  with  constant 
amplitude,  and  a  chirp  with  amplitude  increasing  in  an  RC  time  constant  fashion,  not  to 
exceed  a  given  maximum  value.  These  signals  are  described  below. 

1.  Sinusoidal  Signal 

The  frequency  of  the  sinusoidal  signal  is  selected  randomly  for  every  trial,  where 
the  frequency  range  is  limited  to  avoid  abasing  and  DC  signals.  Figure  5  shows  an 
example  where  f0=  10/512,  and  the  sampling  frequency  fs  =  1 . 


Figure  5.  Sinusoidal  test  signal;  frequency  /„  =10/512,  sampling  frequency  fs  =  1 . 


Constant  Amplitude  Chirp 


The  constant  amplitude  chirp  used  in  the  simulations  is  obtained  by 

f  3470 


s(t)  -  sin 


v  40  +  t  j 


The  signal  sampled  with  sampling  frequency  fs  =  1  is  shown  in  Figure  6.a. 


(4.1) 
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3.  Chirp  With  an  RC  Time  Constant-Like  Amplitude  Increase 

This  signal  has  a  frequency  starting  with  a  high  value  and  decreasing  in  time. 
However,  the  amplitude  of  the  signal  starts  with  a  small  value  and  increases  with  time. 


The  increasing  amplitude  chirp  is  obtained  by 


s(t)  = 


(f +  9)(512-f)  .  (  2;t  547.05  3 
1  -sin 


521 


t  +  35.05 


(4.2) 


This  test  signal  sampled  with  sampling  frequency,  fs  =  1 ,  was  obtained  from  Wavelab 
[13]  and  is  illustrated  in  Figure  6.b. 


B.  SIMULATION  RESULTS 

MATLAB  [14]  simulations  were  performed  to  test  the  performance  of  the 
proposed  schemes.  One  hundred  trials  are  considered  for  both  schemes  for  signal- to- noise 
ratio  (SNR)  values  ranging  between  -6  dB  and  6  dB.  To  perform  the  hypothesis  test  in 
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the  second  step  of  both  FFT-based  and  WT-based  schemes,  the  Bootstrap  MATLAB 
Toolbox  [15]  was  used.  The  MATLAB  code  is  included  in  Appendix  A. 


The  mean  square  error  (MSE)  and  cross- correlation  coefficients  between  original 
and  recovered  signals  were  selected  as  performance  criteria.  The  MSE  was  initially 
selected  as  it  is  commonly  used  in  signal  processing  applications  to  measure  signal 
differences  in  the  time  domain.  It  is  defined  as: 


1  M  f  1  N  |  ^  2 


(4.3) 


where  M  is  the  number  of  trials,  N  is  the  signal  length,  Sj(i)  and  s j(i)  are  the  ith  data 

sample  of  the  original  signal  and  the  recovered  signal  at  f1  trial,  respectively.  Note  that 
the  MSE  may  not  always  be  useful  in  evaluating  actual  performances.  For  example,  in  a 
denoising  scheme  when  the  recovered  signal  is  close  to  the  original  version  in  most  of  the 
signal  duration  but  it  is  very  different  for  a  short  duration,  the  overall  performance  of  the 
scheme  can  be  satisfactory  but  at  the  same  time  the  MSE  may  be  large.  In  addition, 
simulations  showed  the  MSE  performances  to  be  very  similar  on  a  large  portion  of  the 
schemes  investigated.  Therefore,  we  considered  an  additional  performance  criterion, 
based  on  the  cross-correlation  coefficient  to  complement  the  information  given  by  the 
MSE  criterion.  Recall  that  the  normalized  cross- correlation  coefficient  is  commonly  used 
to  evaluate  signal  similarities  and  is  defined  as: 

r  =  - =T - -  £  s(i)* S(i)  ’  (4'4) 

,M2£M  M 

V  1=1  i—  1 


where  ()  denotes  complex  conjugation.  Note  that  the  cross- correlation  coefficient  r 
should  be  equal  to  1  for  a  perfectly  reconstmcted  noise- free  signal.  The  closer  the 
magnitude  of  r  gets  to  0,  the  worse  the  denoised  signal  will  be.  The  distance  measure 
between  1  and  the  normalized  cross- correlation  coefficient  is  given  by 


d 


M 


En-'T. 


7=1 


(4.5) 
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where  M  is  the  number  of  trials  and  rj  is  the  cross- correlation  at  lag  zero  for  the  j,h  trial. 

The  resulting  measure  cl  is  the  additional  criterion  used  to  evaluate  the  denoising 
schemes’  performances  in  extracting  the  noise- tree  signal  from  the  noisy  signal. 

In  our  simulations  sometimes  the  MSE  performances  of  various  schemes  were 
close  to  each  other  and  it  was  difficult  to  make  a  distinction  among  them.  In  this  kind  of 
cases  the  distance  measure  r  was  used  in  making  a  distinction. 

The  proposed  denoising  schemes  are  compared  with  the  original  Wavelet 
shrinkage  algorithm,  which  was  introduced  by  Donoho  and  Johnstone  [16,  17]  to  denoise 
signals  embedded  in  additive  white  Gaussian  noise  with  unit  variance.  Note  that  our 
schemes  do  not  require  knowledge  of  the  noise  variance  levels  to  be  applied,  which  is  not 
the  case  for  the  original  Donoho  and  Johnstone  scheme.  Also  note  that  variable  noise 
variances  had  been  used  in  the  simulations  for  the  proposed  scheme.  Therefore  the  setup 
for  the  simulations  needed  to  be  changed  before  comparing  the  simulations  for  the 
proposed  schemes  with  the  Wavelet  shrinkage  algorithm.  The  Wavelet  shrinkage 
algorithm  is  defined  and  comparison  results  are  presented  in  Appendix  C. 

1.  Fast  Fourier  Transform-Based  Denoising 

Performance  criteria  for  the  three  test  signals  are  illustrated  for  both  thresholding 
options,  in  Figures  7  to  9.  Examples  of  noisy  signals  and  their  recovered  versions 
obtained  by  FT- based  denoising  scheme  are  included  in  Appendix  B. 
a.  Sinusoidal  Signal 

Results,  shown  in  Figure  7,  indicate  that  hard  and  soft  thresholding 
schemes  have  similar  performances  at  all  SNR  values  considered. 
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(a)  Distance  measure  for  sinusoidal  signal 


SNR 

Figure  7.  Distance  measure  d  and  MSE  performance  criteria  for  the  sinusoidal  signal;  FFT-based 

denoising. 


b.  Constant  Amplitude  Chirp 

The  MSE  results  shown  in  Figure  8  indicate  no  significant  differences  in 
thresholding  scheme  performances,  while  the  distance  measure  results  indicate  better 
performance  for  the  soft  thresholding  scheme  at  low  and  medium  SNR  levels. 
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(a)  Distance  measure  for  constant  amplitude  chirp  signal 


SNR 

Figure  8.  Distance  measure  d  and  MSE  performance  criteria  for  constant  amplitude  chirp  signal; 

FFT-based  denoising. 

c.  Chirp  With  an  RC  Time  Constant-Like  Amplitude  Increase 

Results  shown  in  Figure  9  indicate  that  both  thresholding  schemes  have 
similar  MSE  performances  for  SNR  values  below  -3  dB,  and  that  the  soft  thresholding 
scheme  performs  better  for  SNR  values  above  -3  dB.  The  distance  measure  criterion 
indicates  a  better  performance  for  soft  thresholding  for  all  SNR  levels  investigated. 
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(a)  Distance  measure  for  increasing  amplitude  chirp  signal 


SNR 

Figure  9.  Distance  measure  d  and  MSE  performance  criteria  for  increasing  amplitude  chirp  signal; 

FFT-based  denoising. 

2.  Wavelet  Transform-Based  Denoising 

Many  parameters  affect  the  performance  of  the  wavelet  transform-based 
denoising  scheme.  Recall  that  we  consider  four  thresholding  implementation  schemes  on 
the  approximation  coefficients:  approximation  coefficients  left  unperturbed  (method  1), 
approximation  coefficients  thresholded  identical  to  the  detail  coefficients  (method  2), 
approximation  coefficients  median  filtered  with  a  filter  of  length  three  (method  3),  and 
using  an  order- two  predictor  (method  4).  Methods  3  and  4  do  not  perform  as  well  as  the 
first  two  methods  because  they  oversmooth  the  approximation  coefficients  sequence, 
causing  loss  of  signal  power  and  degradation  in  the  performance  of  the  denoising  scheme. 

Other  parameters  are  the  number  of  decomposition  levels  and  the  order  of  the  Daubechies 
wavelet  to  be  used  (3,  4  and  5).  The  numbers  in  the  figure  legends  next  to  the 
thresholding  types  are  the  Daubechies  wavelet  orders.  MSE  values  obtained  for  all 
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schemes  are  close  to  each  other  and  do  not  provide  much  information  about  relative 
performances.  As  a  result,  we  consider  the  distance  measure  d  to  compare  performances. 
In  addition,  we  only  show  results  obtained  for  the  best  of  the  4  methods  for  the 
approximation  coefficients  and  the  best  of  the  3  wavelet  based  decompositions  (3,  4  or  5 
level  decompositions).  Examples  for  noisy  signals  and  their  recovered  versions  obtained 
with  WT- based  denoising  schemes  are  included  in  Appendix  B. 

a.  Sinusoidal  Signal 

Results  show  that  best  performances  are  obtained  with  method  2  for  this 
signal  type.  Figure  10  presents  the  WT-based  scheme  performance  obtained  for  a  4- level 
decomposition  with  different  wavelet  orders  and  both  thresholding  schemes.  Note  that  the 
sinusoidal  signal  has  a  constant  frequency  remaining  in  one  of  the  frequency  bands 
formed  by  the  wavelet  transform  for  the  whole  signal  duration.  Results  show  that  the 
performance  of  the  WT-based  denoising  scheme  depends  on  the  frequency  of  the  specific 
test  signal  and  the  level  of  wavelet  decomposition.  Note  that  test  signal  frequencies  were 
picked  randomly  for  each  trial.  It  should  also  be  noted  that  when  using  method  1,  the 
signal  is  left  untouched  only  when  it  is  located  in  the  lowest  frequency  band  (i.e.,  the 
band  containing  the  approximation  coefficients).  However,  noise  when  present  is  also  left 
untouched  in  that  lowest  frequency  band,  while  no  such  distinction  is  present  in  method 
2. 

Simulation  results  show  that  the  soft  thresholding  scheme  outperforms  the 
hard  thresholding  implementation  at  SNR  values  below  2  dB,  while  all  thresholding 
schemes  perform  similarly  for  higher  SNRs. 
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(a)  Distance  measure  for  sinusoidal  signal,  4  level  decomposition,  method  2 


SNR 

Figure  10.  Distance  measure  and  MSE  for  test  signal  1,  using  a  Wavelet-based  4- level 
decomposition  and  method  2,  with  Daubechies  wavelet  orders  3,  4  and  5. 


b.  Constant  Amplitude  Chirp 

Results  for  the  WT- based  denoising  scheme  using  4- level  decomposition 
and  method  1  are  presented  in  Figure  11.  Method  1  was  shown  here  because  simulations 
indicated  that  method  1  has  significantly  higher  performance  than  method  2  on  the 
constant  amplitude  chirp  signal.  This  result  was  to  be  expected  as  this  signal  has  most  of 
its  power  at  low  frequencies.  Frequency  components  exist  in  several  frequency  bands, 
including  the  one  containing  the  approximation  coefficients.  When  method  2  is  used,  the 
sequence  of  approximation  coefficients  is  tested  for  Gaussianity  and  may  be  estimated  as 
Gaussian  when  the  signal  components  are  very  short,  resulting  in  thresholding  of  the 
signal  components.  However,  such  a  problem  does  not  exist  in  method  1  where  the 
approximation  coefficients  are  left  untouched.  Finally,  simulation  results  show  that  a 
four-level  decomposition  performs  better  than  a  three  or  five-level  decomposition. 
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(a)  Distance  measure  for  constant  amplitude  chirp  signal,  4  level  decomposition,  method  1 


Figure  11.  Performance  for  constant  amplitude  chirp  signal,  using  a  Wavelet  based  4-level 
decomposition  and  method  1,  with  Daubechies  wavelet  orders  3,  4  and  5. 


c.  Chirp  With  an  RC  Time  Constant-Like  Amplitude  Increase 

Simulations  show  method  1  has  the  best  performance  of  the  4  methods 
investigated  and  associated  results  are  presented  in  Figure  12.  Results  shown  in  Figure  12 
indicate  that  similar  performances  are  obtained  for  all  thresholding  schemes  and  wavelet 
orders  considered. 
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(a)  Distance  measure  for  increasing  amplitude  chirp  signal,  4  level  decomposition,  method  1 


-6  -4  -2  0  2  4  6 
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Figure  12.  Performance  for  increasing  amplitude  chirp  signal,  using  a  Wavelet-based  4-level 
decomposition  and  method  1,  with  Daubechies  wavelet  orders  3,  4  and  5. 


In  this  chapter  the  test  signals  were  described  and  the  simulation  results 
were  summarized.  The  conclusions  are  presented  in  the  next  chapter. 
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V.  CONCLUSIONS  AND  FUTURE  WORK 


A.  CONCLUSIONS 

This  work  considered  a  Bootstrap-based  denoising  scheme  to  denoise 
deterministic  signals  embedded  in  additive  Gaussian  noise.  The  Bootstrap  is  a  technique 
derived  for  improving  the  accuracy  of  a  parameter  estimator,  used  especially  in  situations 
where  conventional  techniques  are  not  valid  since  the  data  considered  is  short  in  length. 
This  technique  is  used  for  parameter  estimation  and  hypothesis  testing  in  our  work. 

The  proposed  Bootstrap- based  denoising  scheme  has  four  steps.  First,  the  noisy 
signal  is  transformed  into  a  two-dimensional  domain,  using  the  Short-Time  Fast  Fourier 
or  the  Wavelet  Transform.  Next,  Gaussianity  tests  are  performed  on  the  transformed  data 
on  a  frequency  band  basis.  No  processing  is  applied  to  the  data  when  the  tested  data  is 
found  to  be  non- Gaussian,  as  this  indicates  signal  components  in  the  data  are  dominant 
over  Gaussian  noise  components.  However,  denoising  is  applied  when  the  data  tested  is 
found  to  be  Gaussian,  which  indicates  that  the  noise  components  are  dominant.  Denoising 
is  obtained  by  thresholding  the  frequency  components  in  specific  frequency  bands  to 
minimize  noise  effects.  Finally,  the  recovered  signal  is  obtained  by  applying  the 
appropriate  inverse  transform. 

Denoising  scheme  performances  were  investigated  using  three  test  signals:  a 
sinusoid,  a  constant  amplitude  chirp  and  a  chirp  with  increasing  amplitude.  MSE  and 
cross-correlation  measures  were  selected  to  investigate  the  relative  performances  of  all 
four  denoising  schemes  considered  in  this  work.  In  most  cases,  for  a  given  SNR  the  MSE 
values  obtained  for  schemes  with  different  decomposition  types  and  levels  were  similar  at 
all  SNRs  and  were  not  useful  in  discriminating  between  the  various  schemes  investigated. 
However,  the  distance  measure  was  more  sensitive  and  showed  more  discriminating 
information.  Results  show  that  FFT-based  schemes  perform  better  than  WT-based 
schemes  on  the  stationary  sinusoid  signal  type,  whereas  WT-based  scheme  outperforms 
the  FFT-based  scheme  on  the  chirp  signal  types  considered.  In  the  FFT  case,  soft  and 
hard  thresholding  schemes  perform  similarly  on  the  sinusoid,  while  soft  thresholding 
outperforms  hard  thresholding  on  the  chirp  signals.  Results  also  show  that  soft 
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thresholding  performs  better  on  the  stationary  sinusoidal  signal  type  but  performs  similar 
to  hard  thresholding  on  the  chirps  in  the  WT  case.  Note  that  thresholding  methods  are 
performed  on  the  detail  coefficients  in  the  WT- based  denoising  case.  However,  using 
them  directly  on  the  approximation  coefficients  may  cause  signal  loss.  Therefore,  four 
schemes  were  considered  in  the  simulations:  Leaving  the  approximation  coefficients 
untouched,  using  the  same  thresholding  schemes  as  the  detail  coefficients,  using  a  median 
filter  of  order  3  on  the  approximation  coefficients,  and  using  a  predictor  of  order  2  on  the 
approximation  coefficients.  The  second  method  performs  best  for  the  sinusoidal  test 
signal  whereas  the  first  method  performs  the  best  for  the  chirp  type  test  signals. 

Daubechies  wavelets  of  orders  3,  4  and  5  were  selected  in  the  WT- based 
denoising  scheme.  Simulations  show  that  wavelets  of  this  family  with  higher  orders  did 
not  improve  performance  for  signals  with  the  given  length. 

The  proposed  denoising  scheme  was  compared  with  the  Wavelet  shrinkage 
algorithm,  which  was  introduced  by  Donoho  and  Johnstone  [16,  17]  to  denoise  signals 
embedded  in  additive  white  Gaussian  noise  with  unit  variance.  Results  indicate  that  the 
proposed  WT-based  denoising  scheme  performs  better  than  the  wavelet  shrinkage 
algorithm  for  the  sinusoidal  test  signal.  However,  the  wavelet  shrinkage  algorithm 
outperforms  the  proposed  scheme  for  the  chirp  type  test  signals.  It  should  be  noted  that 
these  results  (Donoho  and  Johnstone)  are  restricted  to  the  case  of  noise  with  unit 
variance. 

B.  FUTURE  WORK 

Our  study  was  restricted  to  one  wavelet  family  (Daubechies)  only  and  further 
investigations  should  consider  other  types  of  wavelets  to  investigate  the  impact  of  the 
specific  wavelet  type  on  the  resulting  denoising  performances. 
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APPENDIX  A.  MATLAB  CODES 


MATLAB  simulations  were  performed  to  were  used  to  test  the  performance  of 
the  denoising  scheme.  The  codes  used  are  presented  in  this  chapter.  The  first  and  second 
sections  present  the  codes  for  the  FFT-  and  WT-based  denoising  schemes,  respectively. 
The  third  section  presents  the  sub -functions  that  were  used  in  both  of  these  codes. 


A.  FFT-BASED  DENOISING 


ooooooooooooooooooooooooooooooooooo 
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o 
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testkurtv7 : 

implement  s 
three  test 

the  FFT-based 
signals 

denoising  scheme 

on  % 

o. 

o 

o. 

SYNTAX  : 

testkurtv7 

o 

o. 

o 

o. 

INPUT  : 

none 

o 

o, 

o 

o. 

OUTPUT  : 

simulation 

results  saved 

on  disk 

o 

o. 

o 

o. 

SUB  FUNC  : 

MakeSignal . m 
kurtosistest.m 

o 

o. 

o 

o. 

o 

o, 

o 

Written  by  Hasan  E.  KAN 


oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


clear 

clc 

close  all 
t=l : 512 ; 

ampl=sqrt ( 10 . A (  [  . 6  0  -.6])); 

kurlim=4 . 5 ; 

for  sigtype=l:3 

for  trial=l:100 

recksoft  =  zeros (length (ampl)  ,  512  )  ; 
reckhard=zeros (length (ampl)  ,  512 )  ; 
if  sigtype==l 

s  =  sin(2*pi*t* (1+ceil (rand*240) ) /512)  ; 
s  =  s/ sqrt (mean(s.A2)  )  ; 

sigtypes= '  sin ' ; 
elseif  sigtype==2 
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s  =  sin  (34. 7. /[. 01:. 01:612]); 
s  =  s (41:552) ; s  =  s / sqrt (mean ( s . "2 )  )  ; 
sigtypes= ' cch ' ; 
elseif  sigtype==3 

s=make signal ( ' Doppler '  ,  521 )  ; 
s  =  s (10:521) ; s  =  s / sqrt (mean ( s . "2 )  )  ; 

sigtypes= ' ech '  ; 

end 

signal=ones (length (ampl )  ,  1 ) *s ; 
sigpwr=sum (signal. ^2, 2) /512; 

n=randn (length (ampl) , 512) ; 

noise= ( (ampl . /std (n ' ) ) ' *ones (1,512)) .  *n; 

noipwr=sum (noise . *2, 2)  /  512; 

SNR=10*logl0 ( sigpwr . / noipwr )  ; 

x=signal+noise; 

for  snr=l : length (ampl) %1 : 21 
f x= zeros (61,32)  ; 
for  segmentno=l : 61 

fx  (segmentno,  : ) =fft (x (snr,  (segmentno- 
1 ) *8  +  1 :  ( segment no- 1 ) *8  +  32)  .*triang(32)  '  )  ; 
end 

fksoft=zeros (61,32)  ; 
f ksof tr=zeros  (61, 17) ; 
f ksof ti=zeros  (61, 17) ; 
fkhard=zeros  (61, 32) ; 
f khardr=zeros  (61, 17) ; 
f khardi=zeros  (61, 17) ; 
kurr= zeros (length (ampl ) , 17 ) ; 
kuri=zeros (length (ampl ) , 17 ) ; 
for  segmentno=l : 17 

[ f ksoftr ( : , segmentno) , f khardr ( : , segmentno) , kurr (segmentno) ] =kurto 
sistest  (real (fx ( : , segmentno) ) , kurlim) ; 

if  imag (fx ( : , segmentno) ) ~=0 

[fksofti ( : , segmentno) , fkhardi ( : , segmentno) , kuri (segmentno) ] =kurto 
sistest ( imag ( f x ( : , segmentno) )  ,  kurlim)  ; 
end 

end 

fksoft ( :, l:17)=fksoftr( :, 1 : 17) + j* fksofti  (  :, 1:17) ; 
fkhard ( : , 1 : 17) =f khardr ( : , 1 : 17) + j*fkhardi ( : ,  1 : 17)  ; 
fksoft (: , 18:32) =con j (fliplr (fksoft ( :,  2:16) ) )  ; 
fkhard ( :, 18:32)=conj (fliplr (fkhard ( : ,  2 : 16)  )  )  ; 
for  segmentno=l : 61 

recksof t ( snr ,  ( segmentno-1 ) *8  +  1 :  (segmentno- 
1 ) *8+32 ) =recksoft ( snr, ( segmentno-1 ) *8+1 : (segmentno- 
1) *8  +  32) +ifft (fksoft (segmentno,  : )  )  ; 
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reckhard ( snr , ( segmentno-1 ) *8+1 : (segmentno- 
1 ) *8+32 ) =reckhard ( snr, ( segmentno-1 ) *8+1 : (segmentno- 
1) *8+32) +if ft ( f khard ( segment no,  : ) )  ; 
end 

end 

recksoft=recksoft*0.5; 
reckhard=reckhard*0 . 5 ; 

msesoft  (:, trial)  =sum  (  ( signal -re  cksoft )  . A2 , 2 ) /512 ; 
msehard  (  : ,  trial )  =sum  (  ( signal -reckhard)  . A2 , 2 ) /512 ; 
for  snr=l : 2 1 

corsoft (snr, trial) =xcorr (s, recksoft (snr, : ) , 0, ' coef f ' ) ; 

cor hard (snr, trial) =xcorr (s, reckhard (snr,  :)  ,0,  ' coef f '  )  ; 
end 

end 

save ( sprint f ( ' mse%s ' , sigtypes ) ,  ' msemed '  ,  'msesoft',  ' msehard ' ) 
save ( sprint f ( ' cor%s ' , sigtypes ) ,  'corsoft',  ' corhard '  ) 

end 


B.  WT-BASED  DENOISING 
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testkurtv7 : 

implements  the  WT -based  denoising  scheme  on 

o, 

o 

three  test  signals 

o. 

o 

o. 

SYNTAX  : 

testwav2 

o 

o. 

o 

o. 

INPUT  : 

none 

o 

o. 

o 

o. 

OUTPUT  : 

simulation  results  saved  on  disk 

o 

o, 

o 

o. 

SUB  FUNC  : 

MakeSignal . m 

o 

o. 

o 

kurtosistest .m 

o, 

o 

o. 

o 

Written  by  Hasan  E.  KAN  % 

o. 

o 


oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 


clear 

clc 

close  all 

t=l : 512 ; 

ampl=sqrt ( 10 . A (  [  . 6  0  -.6])); 
kurlim=4 . 5 ; 

rec5db3soft=zeros ( length (ampl ) , 512 , 4 ) ; 
rec5db3hard= zeros (length (ampl )  ,  512 , 4  )  ; 
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rec4db3soft=zeros (length (ampl)  ,512,4)  ; 
rec4db3hard= zeros (length (ampl)  ,512,4)  ; 
rec3db3soft  =  zeros (length (ampl)  ,512,4); 
rec3db3hard=zeros (length (ampl)  ,512,4)  ; 

rec5db4soft  =  zeros ( length (ampl )  ,  512 ,  4  )  ; 
rec5db4 ha rd= zeros (length (ampl)  ,512,4); 
rec4db4soft  =  zeros (length (ampl )  ,  512 ,  4  )  ; 
rec4db4 ha rd= zeros (length (ampl)  ,512,4)  ; 
rec3db4soft  =  zeros ( length (ampl )  ,  512 ,  4 )  ; 
rec3db4hard=zeros (length (ampl)  ,512,4)  ; 

rec5db5soft  =  zeros (length (ampl)  ,512,4); 
rec5db5 ha rd= zeros (length (ampl)  ,512,4)  ; 
rec4db5soft  =  zeros (length (ampl)  ,512,4)  ; 
rec4db5 ha rd= zeros (length (ampl)  ,512,4); 
rec3db5soft  =  zeros (length (ampl)  ,512,4); 
rec3db5hard=zeros (length (ampl)  ,512,4); 

for  sigtype=l:3 

for  trial=l:100 
if  sigtype==l 

s  =  sin ( 2  *pi*t * (1+ceil ( rand* 2 4 0 ))/512);s  =  s/ sqrt (sum (s . A2 ) / 5 1 2 )  ; 
sigtypes= ' sin ' ; 
elseif  sigtype==2 

s=sin(34.7./ [ .01: .01:612] ) ;s=s (41:552) ; s=s/ sqrt ( sum (s.A2) /512) ; 
sigtypes= ' cch ' ; 
elseif  sigtype==3 

s=make signal ( ' Doppler ',521) ;s  =  s (10:521) ; s=s/ sqrt  ( sum  (s.A2)/512)  ; 
sigtypes= ' ech ' ; 

end 

signal=ones (length (ampl )  ,  1 ) *s ; 
sigpwr=sum (signal. A2, 2) /512; 

n=randn (length (ampl) , 512) ; 

noise= ( (ampl . /std (n ' ) ) ' *ones (1,512)) . *n; 

noipwr=sum (noise. A2, 2) /512; 

SNR=10*logl0 ( sigpwr . / noipwr )  ; 

x=signal+noise; 

for  snr=l : length (ampl) 

[cdb3 (snr,  :),ldb3] =wavedec (x (snr,  : )  ,  3,  ' db3 ' ) ; 
a3db3 (snr, : ) =appcoef (cdb3 (snr, : ) , ldb3, ' db3 ' , 3) ; 

[dldb3 ( snr , : ) , d2db3 ( snr , : ) , d3db3 ( snr , : ) ] =detcoef ( cdb3 ( snr , : ) , ldb3 
, [1,2,3] ) ; 

[a4db3(snr, :),d4db3(snr, :)] =dwt (a3db3(snr, :), ' db3 ' ) ; 
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[a5db3 (snr, :) ,d5db3 (snr, :) ] =dwt (a4db3 (snr, :  )  , ' db3 ' ) ; 


[ cdb4 (snr, : ) , ldb4] =wavedec ( x ( snr , :),3,'db4'); 
a3db4 (snr,  : ) =appcoef (cdb4 (snr,  :),ldb4,  ' db4 ' , 3) ; 

[dldb4 (snr, : ) , d2db4 (snr, : ) , d3db4 (snr, : ) ] =detcoef ( cdb4 ( snr , : ) , ldb4 
, [1,2,3] ) ; 

[a4db4 (snr, :) ,d4db4 (snr, :) ] =dwt (a3db4 (snr, :) , ' db4 ' ) ; 
[a5db4 (snr, :) ,d5db4 (snr, :) ] =dwt (a4db4(snr, :), ' db4 ' ) ; 

[ cdb5 (snr, : ) , ldb5] =wavedec ( x ( snr , : ) , 3 , ' db5 ' ) ; 
a3db5 (snr,  : ) =appcoef (cdb5 (snr,  :),ldb5,  ' db5 '  ,  3)  ; 

[dldb5 ( snr , : ) , d2db5 ( snr , : ) , d3db5 ( snr , : ) ] =detcoef ( cdb5 (snr, : ) , ldb5 
, [1,2,3] ) ; 

[a4db5(snr, :),d4db5(snr, :)] =dwt (a3db5 (snr, : ) , ' db5 ' ) ; 
[a5db5(snr, :) ,d5db5 (snr, :) ] =dwt (a4db5(snr, :), ' db5 ' ) ; 

end 

ra5db3sof t=zeros (length (ampl) , length (a5db3) ) ; ra5db3hard=zeros (len 
gth (ampl) , length (a5db3) )  ; 

ra4db3sof t=zeros (length (ampl) , length (a4db3) ) ; ra4db3hard=zeros (len 
gth (ampl) , length (a4db3)  )  ; 

ra3db3sof t=zeros (length (ampl) , length (a3db3) ) ; ra3db3hard=zeros (len 
gth (ampl) , length (a3db3) )  ; 

rd5db3soft=zeros (length (ampl) , length (d5db3 ( 1 , :))); rd5db3hard=zero 
s  ( length (ampl ) , length (d5db3 ( 1 ,  :  )  )  )  ; 

rd4db3sof t  =  zeros (length (ampl)  ,  length (d4db3) )  ;  rd4db3hard=zeros (len 
gth (ampl) , length (d4db3)  )  ; 

rd3db3sof t=zeros (length (ampl) , length (d3db3) ) ; rd3db3hard=zeros (len 
gth (ampl) , length (d3db3) ) ; 

rd2db3sof t=zeros (length (ampl)  ,  length (d2db3) )  ;  rd2db3hard=zeros (len 
gth (ampl) , length (d2db3)  )  ; 

rdldb3sof t=zeros (length (ampl) , length (dldb3) ) ; rdldb3hard=zeros (len 
gth (ampl) , length (dldb3) ) ; 


ra5db4soft=zeros (length (ampl ) , length (a5db 4) ) ; ra5db4hard= zeros (len 
gth (ampl) , length (a5db4 ) ) ; 

ra4db4soft=zeros (length (ampl ) , length (a4db 4) ) ; ra4db4hard=zeros (len 
gth (ampl) , length (a4db4 )  )  ; 

ra3db4soft=zeros (length (ampl ) , length (a3db 4) ) ; ra3db4hard= zeros (len 
gth (ampl) , length (a3db4 ) ) ; 
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rd5db4soft=zeros (length (ampl ) , length (d5db 4) ) ; rd5db4hard=zeros ( len 
gth (ampl) , length (d5db4 ) ) ; 

rd4db4soft=zeros (length (ampl ) , length (d4db 4) ) ; rd4db4hard=zeros (len 
gth (ampl) , length (d4db4 )  )  ; 

rd3db4soft=zeros (length (ampl ) , length (d3db 4) ) ; rd3db4hard=zeros (len 
gth (ampl) , length (d3db4 )  )  ; 

rd2db4soft=zeros (length (ampl ) , length (d2db 4 ) ) ; rd2db4hard=zeros (len 
gth (ampl) , length (d2db4 ) ) ; 

rdldb4soft  =  zeros (length (ampl )  ,  length (dldb 4) )  ;  rdldb4hard=zeros (len 
gth (ampl) , length (dldb 4 )  )  ; 

ra5db5sof t=zeros (length (ampl) , length (a5db5) ) ; ra5db5hard=zeros (len 
gth (ampl) , length (a5db5)  )  ; 

ra4db5sof t=zeros (length (ampl) , length (a4db5) ) ; ra4db5hard=zeros (len 
gth (ampl) , length (a4db5) ) ; 

ra3db5sof t=zeros (length (ampl) , length (a3db5) )  ;  ra3db5hard=zeros (len 
gth (ampl) , length (a3db5) )  ; 

rd5db5sof t=zeros (length (ampl) , length (d5db5) ) ; rd5db5hard=zeros (len 
gth (ampl) , length (d5db5) ) ; 

rd4db5sof t  =  zeros (length (ampl)  ,  length (d4db5) )  ;  rd4db5hard=zeros (len 
gth (ampl) , length (d4db5)  )  ; 

rd3db5sof t=zeros (length (ampl) , length (d3db5) ) ; rd3db5hard=zeros (len 
gth (ampl) , length (d3db5) ) ; 

rd2db5sof t=zeros (length (ampl) , length (d2db5) ) ; rd2db5hard=zeros (len 
gth (ampl) , length (d2db5) ) ; 

rdldb5sof t=zeros (length (ampl) , length (dldb5) ) ; rdldb5hard=zeros (len 
gth (ampl) , length (dldb5)  )  ; 
for  method=l : 4 

if  method==l 

ra5db3sof t=a5db3 ; ra5db3hard=a5db3 ; 
ra4db3sof t=a4db3 ; ra4db3hard=a4db3 ; 
ra3db3sof t=a3db3 ; ra3db3hard=a3db3 ; 
ra5db4  sof t=a5db4 ; ra5db4hard=a5db4 ; 
ra4db4  sof t=a4db4 ; ra4db4hard=a4db4 ; 
ra3db4  sof t=a3db4 ; ra3db4hard=a3db4 ; 
ra5db5sof t=a5db5; ra5db5hard=a5db5 ; 
ra4db5sof t=a4db5; ra4db5hard=a4db5; 
ra3db5sof t=a3db5 ; ra3db5hard=a3db5 ; 
methodl='no  operation  on  approximation  coeff.'; 
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elseif  method==2 

for  snr=l : length (ampl ) 

[ra5db3soft (snr,  : ) , ra5db3hard (snr,  : ) , kura5db3 (snr) ] =kurtosistest  ( 
a5db3 (snr, : ) , kurlim) ; 

[ra4db3soft (snr, : ) , ra4db3hard ( snr , : ) , kura4db3 (snr) ] =kurtosistest ( 
a4db3 (snr,  : ) , kurlim)  ; 

[ra3db3 soft (snr, : ) , ra3db3hard ( snr , : ) , kura3db3 ( snr ) ] =kurtosistest ( 
a3db3 (snr,  : ) , kurlim)  ; 

[ra5db4 soft (snr,  : ) , ra5db4hard ( snr ,  : ) , kura5db4 (snr ) ] =kurtosistest  ( 
a5db4 (snr, : ) , kurlim) ; 

[ra4db4soft (snr, : ) , ra4db4hard (snr, :) ,kura4db4 (snr) ] =kurtosistest ( 
a4db4 (snr,  :) , kurlim)  ; 

[ra3db4soft  (snr,  : ) , ra3db4hard ( snr ,  : ) , kura3db4 (snr) ] =kurtosistest ( 
a3db4 (snr,  : ) , kurlim)  ; 

[ra5db5soft (snr,  : ) , ra5db5hard (snr,  : ) , kura5db5 (snr) ] =kurtosistest  ( 
a5db5 (snr, : ) , kurlim) ; 

[ra4db5soft (snr, : ) , ra4db5hard (snr, : ) , kura4db5 (snr) ] =kurtosistest ( 
a4db5 (snr,  : ) , kurlim)  ; 

[ra3db5soft  (snr,  : ) , ra3db5hard (snr,  : ) , kura3db5 (snr) ] =kurtosistest ( 
a3db5 (snr,  : ) , kurlim)  ; 

end 

methodl= ' using  threshold  on  approximation 

coef f .  '  ; 

elseif  method==3 

ra5db3sof t=medf ilt 1 (a5db3 ' , 3 ) ra5db3hard=ra5db3sof t ; 

ra4db3sof t=medf ilt 1 (a4db3 ' , 3 )  ' ; ra4db3hard=ra4db3sof t ; 

ra3db3sof t=medf ilt 1 (a3db3 ' , 3 ) ' ; ra3db3hard=ra3db3sof t ; 

ra5db4 sof t=medf iltl (a5db4 ' , 3) ' ; ra5db4hard=ra5db4 sof t ; 

ra4db4 sof t=medf iltl (a4db4 ' , 3) ' ; ra4db4hard=ra4db4 sof t ; 

ra3db4 sof t=medf iltl (a3db4 ' , 3)  ' ; ra3db4hard=ra3db4 sof t ; 

ra5db5sof t=medf ilt 1 (a5db5 ' , 3 ) ra5db5hard=ra5db5sof t ; 

ra4db5sof t=medf ilt 1 (a4db5 ' , 3 ) ' ; ra4db5hard=ra4db5sof t ; 

ra3db5sof t=medf ilt 1 (a3db5 ' , 3 )  ra3db5hard=ra3db5sof t ; 

methodl= ' using  medfilt3  on  approximation  coeff.'; 
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else 


rra5db3=xcorr (a5db3 ' , 2 , ' unbiased ' ) ;rra5db3=rra5db3 (3:5, 4  * [0:2] +1) 

r 

rra4db3=xcorr (a4db3',2,  'unbiased') ;rra4db3=rra4db 3(3:5,  4* [ 0 : 2 ] +1 ) 

r 

rra3db3=xcorr (a3db3 ' , 2 , ' unbiased ' ) ; rra3db3=rra3db3 (3:5, 4* [0:2] +1) 

r 

rra5db4=xcorr (a5db4',2, 'unbiased') ;rra5db4=rra5db4 (3:5, 4* [0:2] +1) 

r 

rra4db4=xcorr (a4db4',2, 'unbiased') ;rra4db4=rra4db4 (3:5, 4* [0:2] +1) 

r 

rra3db4=xcorr (a3db4',2, 'unbiased') ;rra3db4=rra3db4 (3:5, 4* [0:2] +1) 

r 

rra5db5=xcorr (a5db5 ' , 2 , ' unbiased ' ) ; rra5db5=rra5db5 (3:5, 4* [0:2] +1) 

r 

rra4db5=xcorr (a4db5',2,  'unbiased') ;rra4db5=rra4db 5(3:5, 4  * [ 0 : 2 ] +1 ) 

r 

rra3db5=xcorr (a3db5 ' , 2 , ' unbiased ' ) ; rra3db5=rra3db5 (3:5, 4* [0:2] +1) 

r 

for  snr=l : length (ampl) 

Rxa5db3=toeplit z(rra5db3(:,snr)  ) ; Rxxa5db3=Rxa5db3 ( 1 : 2 ,  1 : 2 )  ; 

Rxa4db3=toeplit z (rra4db3 ( : , snr )  ) ; Rxxa4db3=Rxa4db3 ( 1 : 2 ,  1 : 2 )  ; 

Rxa3db3=toeplit  z (rra3db3 ( : , snr) ) ; Rxxa3db3=Rxa3db3 ( 1 : 2 , 1 : 2 ) ; 

Rxa5db4=toeplit  z (rra5db4 ( : , snr)  ) ; Rxxa5db4=Rxa5db4 ( 1 : 2 , 1 : 2 ) ; 

Rxa4db4=toeplit  z (rra4db4 ( : , snr )  ) ; Rxxa4db4=Rxa4db4 ( 1 : 2 ,  1 : 2 ) ; 

Rxa3db4=toeplit z (rra3db4 ( : , snr)  ) ; Rxxa3db4=Rxa3db4 ( 1 : 2 ,  1 : 2 )  ; 

Rxa5db5=toeplitz (rra5db5 ( : , snr)  ) ; Rxxa5db5=Rxa5db5 (1 : 2,  1 : 2)  ; 

Rxa4db5=toeplit  z (rra4db5 ( : , snr )  ) ; Rxxa4db5=Rxa4db5 ( 1 : 2 , 1 : 2 ) ; 

Rxa3db5=toeplit  z (rra3db5 ( : , snr ) ) ; Rxxa3db5=Rxa3db5 ( 1 : 2 ,  1 : 2 ) ; 

aa5db3=Rxxa5db3\Rxa5db3 (2:3,1) ; 
aa4db3=Rxxa4db3\Rxa4db3 (2:3,1) ; 
aa3db3=Rxxa3db3\Rxa3db3 (2:3,1) ; 
aa5db4=Rxxa5db4 \Rxa5db4 (2:3,1) ; 
aa4db4=Rxxa4db4\Rxa4db4 (2:3,1) ; 
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aa3db4=Rxxa3db4 \Rxa3db4 (2:3,1)  ; 
aa5db5=Rxxa5db5\Rxa5db5 (2:3,1)  ; 
aa4db5=Rxxa4db5\Rxa4db5 (2:3,1)  ; 
aa3db5=Rxxa3db5\Rxa3db5 (2:3,1)  ; 
ra5db3soft(snr,  :)=filter([0; 
aa5db3] , 1 , a5db3 (snr,  :)  )  ; 

ra4db3soft(snr,  : ) =f ilter ( [ 0 ; 
aa4db3] , 1 , a4db3 (snr,  :)  )  ; 

ra3db3soft(snr,  :)=filter([0; 
aa3db3] , 1 , a3db3 (snr,  :)  )  ; 

ra5db4soft(snr,  :)=filter([0; 
aa5db4] , 1 , a5db4 (snr,  :)  )  ; 

ra4db4soft(snr,  : ) =f ilter ( [ 0 ; 
aa4db4] , 1 , a4db4 (snr,  :)  )  ; 

ra3db4soft(snr,  : ) =f ilter ( [ 0 ; 
aa3db4] , 1 , a3db4 (snr,  :)  )  ; 

ra5db5soft(snr,  :)=filter([0; 
aa5db5] , 1 , a5db5 (snr,  :)  )  ; 

ra4db5soft(snr,  :)=filter([0; 
aa4db5],l,a4db5(snr,  :))  ; 

ra3db5soft(snr,  :)=filter([0; 
aa3db5] , 1 , a3db5 (snr,  :)  )  ; 

end 

ra5db3soft(:,l:2)=a5db3(:,l:2) ; 
ra4db3soft ( :, I:2)=a4db3  (  : , 1:2) ; 
ra3db3soft ( :, I:2)=a3db3  (  : , 1:2) ; 
ra5db4soft ( :, I:2)=a5db4 ( :, 1:2) ; 
ra4db4 sof t ( : , 1 : 2 ) =a4db4 ( : ,  1 : 2 )  ; 
ra3db4 sof t ( : , 1 : 2 ) =a3db4 ( : , 1 : 2 )  ; 
ra5db5soft(:,l:2)=a5db5(:,l:2)  ; 
ra4db5soft ( :, I:2)=a4db5 ( :,  1:2)  ; 
ra3db5soft ( :, I:2)=a3db5 ( :,  1:2)  ; 
ra5db3hard=ra5db3sof t ; 
ra4db3hard=ra4db3sof t ; 
ra3db3hard=ra3db3sof t ; 
ra5db4hard=ra5db4  sof t ; 
ra4db4hard=ra4db4  sof t ; 
ra3db4hard=ra3db4  sof t ; 
ra5db5hard=ra5db5sof t ; 
ra4db5hard=ra4db5sof t ; 
ra3db5hard=ra3db5sof t ; 

methodl= ' prediction  on  approximation  coeff.'; 
end 

for  snr=l : length (ampl ) 

[rd5db3soft (snr, : ) , rd5db3hard (snr, : ) , kurd5db3 (snr) ] =kurtosistest ( 
d5db3 (snr,  : ) , kurlim)  ; 

[rd4db3soft (snr, : ) , rd4db3hard (snr, : ) ,kurd4db3 (snr) ] =kurtosistest ( 
d4db3 (snr,  :) , kurlim)  ; 
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[rd3db3soft 

(snr, : ) , rd3db3hard i 

(snr, 

:)  , 

kurd3db3 

(snr) ] 

=kurtosistest 

d3db3 ( snr , : 

) , kurlim) ; 

[rd2db3soft 

(snr ,  : )  ,  rd2db3hard i 

(snr, 

:)  , 

kurd2db3 

(snr) ] 

=kurtosistest 

d2db3 (snr, : 

) , kurlim) ; 

[rdldb3soft 

(snr, : ) , rdldb3hard i 

(snr, 

:)  , 

kurdldb3 

(snr) ] 

=kurtosistest 

dldb3 (snr, : 

) , kurlim) ; 

[ rd5db4  sof t 

( snr ,  : )  ,  rd5db4hard i 

(snr, 

:)  , 

kurd5db4 

(snr) ] 

=kurtosistest 

d5db4 (snr, : 

) , kurlim) ; 

[ rd4db4  sof t 

( snr , : ) , rd4db4hard i 

(snr, 

:)  , 

kurd4db4 

(snr) ] 

=kurtosistest 

d4db4 (snr, : 

) , kurlim) ; 

[ rd3db4  sof t 

(snr, : ) , rd3db4hard i 

(snr. 

:)  , 

kurd3db4 

(snr) ] 

=kurtosistest 

d3db4 (snr, : 

) , kurlim) ; 

[ rd2db4  sof t 

( snr ,  : )  ,  rd2db4hard i 

(snr, 

:)  , 

kurd2db4 

(snr) ] 

=kurtosistest 

d2db4 (snr, : 

) , kurlim) ; 

[ rdldb4  sof t 

( snr , : ) , rdldb4hard i 

(snr, 

:)  , 

kurdldb4 

(snr) ] 

=kurtosistest 

dldb4 (snr, : 

) , kurlim) ; 

[rd5db5soft 

(snr, : ) , rd5db5hard i 

(snr, 

:)  , 

kurd5db5 

(snr) ] 

=kurtosistest 

d5db5 (snr, : 

) , kurlim) ; 

[rd4db5soft 

( snr ,  : )  ,  rd4db5hard i 

(snr, 

:)  , 

kurd4db5 

(snr) ] 

=kurtosistest 

d4db5 (snr, : 

) , kurlim) ; 

[rd3db5soft 

( snr , : ) , rd3db5hard i 

(snr, 

:)  , 

kurd3db5 

(snr) ] 

=kurtosistest 

d3db5 (snr, : 

) , kurlim) ; 

[rd2db5soft 

(snr, : ) , rd2db5hard i 

(snr, 

:)  , 

kurd2db5 

(snr) ] 

=kurtosistest 

d2db5 (snr, : 

) , kurlim) ; 

[rdldb5soft 

(snr, : ) , rdldb5hard i 

(snr, 

:)  , 

kurdldb5 

(snr) ] 

=kurtosistest 

dldb5 (snr, : 

) , kurlim) ; 

end 

R5db3soft= [ra5db3soft 

rd5db3sof t 

rd4db3sof t 

rd3db3sof t 

rd2db3soft  rdldb3soft]; 

R5db3hard= [ra5db3hard 

rd5db3hard 

rd4db3hard 

rd3db3hard 

rd2db3hard  rdldb3hard] ; 

R4db3soft= [ra4db3soft 

rd4db3sof t 

rd3db3sof t 

rd2db3sof t 

rdldb3sof t ] ; 

R4db3hard= [ra4db3hard 

rd4db3hard 

rd3db3hard 

rd2db3hard 

rdldb3hard] ; 

R3db3soft= [ra3db3soft 

rd3db3soft 

rd2db3sof t 

rdldb3sof t 
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rdldb3hard] 

R3db3hard= [ra3db3hard 

r 

rd3db3hard 

rd2db3hard 

rd3db4  sof t 

R5db4soft=[ ra5db4  sof t 
rd2db4soft  rdldb4soft]; 

rd5db4  sof t 

rd4db4  sof t 

rd3db4hard 

R5db4hard= [ra5db4hard 
rd2db4hard  rdldb4hard] ; 

rd5db4hard 

rd4db4hard 

rd2db4  sof t 

R4db4soft=[ ra4db4  sof t 
rdldb4  sof t ] ; 

rd4db4  sof t 

rd3db4  sof t 

rd2db4hard 

R4db4hard= [ra4db4hard 
rdldb4hard] ; 

rd4db4hard 

rd3db4hard 

rdldb4  sof t ] 

R3db4soft=[ ra3db4  sof t 

r 

rd3db4  sof t 

rd2db4  sof t 

rdldb4hard] 

R3db4hard= [ra3db4hard 

r 

rd3db4hard 

rd2db4hard 

rd3db5sof t 

R5db5soft= [ra5db5soft 
rd2db5soft  rdldb5soft]; 

rd5db5sof t 

rd4db5sof t 

rd3db5hard 

R5db5hard= [ra5db5hard 
rd2db5hard  rdldb5hard] ; 

rd5db5hard 

rd4db5hard 

rd2db5sof t 

R4db5soft= [ra4db5soft 
rdldb5sof t ] ; 

rd4db5sof t 

rd3db5sof t 

rd2db5hard 

R4db5hard= [ra4db5hard 
rdldb5hard] ; 

rd4db5hard 

rd3db5hard 

rdldb5sof t ] 

R3db5soft= [ra3db5soft 

r 

rd3db5sof t 

rd2db5sof t 

R3db5hard= [ra3db5hard 

rd3db5hard 

rd2db5hard 

rdldb5hard] ; 

for  snr=l : length (ampl) 

rec5db3sof t ( snr,  method) =waverec (R5db3soft ( snr,  :),[ 20  20  36  68 
131  258  512] , ' db3 ' ) ; 

rec5db3hard (snr,  : , method) =waverec (R5db3hard (snr,:),  [20  20 

36  68  131  258  512] , 'db3 ' ) ; 

rec4db3sof t ( snr ,:, method) =waverec (R4db3soft ( snr ,:),[ 36  36  68  131 
258  512] , ' db3 ' ) ; 

rec4db3hard (snr, : , method) =waverec (R4db3hard (snr,:), [36  36 

68  131  258  512] , 'db3 ' ) ; 

rec3db3soft ( snr,  :, method) =waverec (R3db3soft ( snr,  :),[  68  68  131  258 
512] , ' db3 ' ) ; 

rec3db3hard (snr, : , method) =waverec (R3db3hard (snr,:), [68  68 

131  258  512] , ' db3 ' ) ; 


rec5db4 soft ( snr ,:, method) =waverec (R5db4 soft ( snr ,:),[ 22  22  38  70 
133  259  512] , 'db4 ' ) ; 
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rec5db4hard ( snr , :, method) =waverec (R5db4hard ( snr , :), [22  22 

38  70  133  259  512] ,  'db4 ' )  ; 

rec4db4 soft ( snr method) =waverec (R4db4 soft  ( snr 38  38  70  133 
259  512] , ' db4 ' ) ; 

rec4db4hard ( snr , method) =waverec (R4db4hard ( snr , :), [38  38 

70  133  259  512] ,  'db4 ' )  ; 

rec3db4 soft ( snr method) =waverec (R3db4 soft ( snr ,:),[ 70  70  133  259 
512] , ' db4 ' ) ; 

rec3db4hard ( snr ,  : , method) =waverec (R3db4hard (snr,:),  [70  70 

133  259  512]  ,  'db4 '  )  ; 


rec5db5sof t ( snr ,:, method) =waverec (R5db5soft ( snr ,:),[ 24  24  40  71 
134  260  512]  ,  ' db5 ' ) ; 

rec5db5hard ( snr ,  :, method) =waverec (R5db5hard ( snr ,  :),  [24  24 

40  71  134  260  512]  ,  'db5 ' ) ; 

rec4db5sof t ( snr ,:, method) =waverec (R4db5soft ( snr ,:),[ 4 0  40  71  134 
260  512] ,  ' db5 '  )  ; 

rec4db5hard ( snr ,  :, method) =waverec (R4db5hard ( snr ,  :),  [40  40 

71  134  260  512]  ,  'db5 ' ) ; 

rec3db5soft ( snr ,:, method) =waverec (R3db5soft ( snr ,:),[ 7 1  71  134  260 
512] , ' db5 ' ) ; 

rec3db5hard ( snr ,  :, method) =waverec (R3db5hard ( snr ,  :),  [71  71 

134  260  512]  ,  ' db5 '  )  ; 

end 


end 

signa=repmat (signal ,  [1,1,4]); 


ms e5db3soft(:, trial, :) =sum ( 
rec5db3soft) .A2,2)/512; 

mse5db3hard ( : , trial, : ) =sum ( 
rec5db3hard) . A2,2)/512; 

ms e4db3soft(:, trial, :) =sum ( 
rec4db3soft) .A2,2)/512; 

mse4db3hard (:, trial, :)=sum( 
rec4db3hard) .A2,2)/512; 

ms e3db3soft(:, trial, :) =sum ( 
rec3db3soft) .A2,2)/512; 

mse3db3hard ( : , trial, : ) =sum ( 
rec3db3hard) . A2,2)/512; 


( signa- 
( signa- 
( signa- 
( signa- 
( signa- 
( signa- 
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ms e5db4soft(:, trial, :) =sum ( ( signa- 
rec5db4soft) .A2,2)/512; 

mse5db4hard (:, trial, : ) =sum ( (signa- 
rec5db4hard) . A2,2)/512; 

ms e4db4soft(:, trial, :) =sum ( ( signa- 
rec4db4soft) .A2,2)/512; 

mse4db4hard (:, trial, : ) =sum ( (signa- 
rec4db4hard) .A2,2)/512; 

ms e3db4soft(:, trial, :) =sum ( ( signa- 
rec3db4soft) .A2,2)/512; 

mse3db4hard (:, trial, : ) =sum ( (signa- 
rec3db4hard) .A2,2)/512; 

ms e5db5soft(:, trial, :) =sum ( (signa- 
rec5db5soft) .A2,2)/512; 

mse5db5hard (:, trial, : ) =sum ( ( signa- 
rec5db5hard) .A2,2)/512; 

ms e4db5soft(:, trial, :) =sum ( (signa- 
rec4db5soft) .A2,2)/512; 

mse4db5hard (:, trial, : ) =sum ( ( signa- 
rec4db5hard) .A2,2)/512; 

ms e3db5soft(:, trial, :) =sum ( (signa- 
rec3db5soft) .A2,2)/512; 

mse3db5hard (:, trial, : ) =sum ( (signa- 
rec3db5hard) .A2,2)/512; 
for  snr=l : 21 

for  method=l : 4 

cor3db3soft (snr,  trial, method) =xcorr(s,rec3db3soft(snr,  :, method) , 0 
, ' coef f ' ) ; 

cor3db4soft (snr, trial, method) =xcorr (s,rec3db4soft(snr, :, method) , 0 
, ' coef f ' ) ; 

cor3db5soft (snr, trial, method) =xcorr(s,rec3db5soft(snr,  :,  method)  ,  0 
, ' coef f ' ) ; 

cor4db3soft (snr, trial, method) =xcorr (s,rec4db3soft(snr, :, method) , 0 
, ' coef f ' ) ; 

cor4db4soft (snr, trial, method) =xcorr (s,rec4db4soft(snr,  :, method) , 0 
, ' coef f ' ) ; 

cor4db5soft (snr, trial, method) =xcorr (s,rec4db5soft(snr, :, method) , 0 
, ' coef f ' ) ; 

cor5db3soft (snr, trial, method) =xcorr(s,rec5db3soft(snr, :, method) , 0 
, ' coef f ' ) ; 

cor5db4soft (snr, trial, method) =xcorr (s,rec5db4soft(snr,  :, method)  ,  0 
, ' coef f ' ) ; 
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cor5db5soft (snr, trial, method) =xcorr(s,rec5db5soft(snr,  :  ,  method)  ,  0 
, ' coef f ' ) ; 

cor3db3hard (snr, trial, method) =xcorr ( s , rec3db3hard (snr, : , method) , 0 
, ' coef f ' ) ; 

cor3db4hard (snr, trial, method) =xcorr ( s ,  rec3db4hard (snr,  : , method) , 0 
, ' coef f ' ) ; 

cor3db5hard (snr, trial, method) =xcorr ( s , rec3db5hard (snr, : , method) , 0 
, ' coef f ' ) ; 

cor4db3hard (snr, trial, method) =xcorr ( s , rec4db3hard (snr, : , method) , 0 
, ' coef f ' ) ; 

cor4db4hard (snr, trial, method) =xcorr ( s , rec4db4hard (snr,  :  ,  method)  ,  0 
, ' coef f ' ) ; 

cor4db5hard (snr, trial, method) =xcorr ( s , rec4db5hard (snr, : , method) , 0 
, ' coef f ' ) ; 

cor5db3hard (snr, trial, method) =xcorr ( s ,  rec5db3hard (snr,  : , method) , 0 
, ' coef f ' ) ; 

cor5db4hard (snr, trial, method) =xcorr ( s , rec5db4hard (snr, : , method) , 0 
, ' coef f ' ) ; 

cor5db5hard (snr, trial, method) =xcorr ( s , rec5db5hard (snr, : , method) , 0 
, ' coef f ' ) ; 

end 

end 

end 


save (sprintf ( 'mse%s31vl ' , sigtypes) 
3hard', ' mse3db4 sof t ' , ' mse3db4hard ' 

save (sprintf ( 'mse%s41vl ' , sigtypes) 
3hard', ' mse4db4 sof t ' , ' mse4db4hard ' 


'msemed' , 'mse3db3soft ' , 'mse3db 
'mse3db5soft ' , ' mse3db5hard ' ) 

'msemed', ' mse4db3sof t ' , 'mse4db 
'mse4db5soft ' , ' mse4db5hard ' ) 


save (sprintf ( 'mse%s51vl ' , sigtypes) , 'msemed' , 'mse5db3soft ' , 'mse5db 
3hard', ' mse5db4 sof t ' , ' mse5db4hard ' , ' mse5db5soft ' , ' mse5db5hard ' ) 

save (sprintf ( ' cor%s31vl ', sigtypes ) , ' cor 3db3 soft ' , ' cor3db3hard ' , ' c 
or3db4soft ' , ' cor3db4hard' , ' cor3db5soft ' , ' cor3db5hard' ) 


save (sprintf ( ' cor%s41vl ' , sigtypes) , ' cor4db3soft ' , ' cor4db3hard' , ' c 
or4db4soft ' , ' cor4db4hard' , ' cor4db5soft ' , ' cor4db5hard' ) 


save (sprintf ( ' cor%s51vl ', sigtypes), 'cor5db3soft', ' cor5db3hard ' , ' c 

or5db4soft ' , ' cor5db4hard' , ' cor5db5soft ' , ' cor5db5hard' ) 

end 
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C.  SUB  FUNCTIONS 
1.  MakeSignal 

This  code  is  part  of  Wavelab  version  802  [11]. 


function  sig  =  MakeSignal (Name, n) 

%  MakeSignal  --  Make  artificial  signal 
%  Usage 

%  sig  =  MakeSignal (Name, n) 

%  Inputs 

%  Name  string:  'HeaviSine',  'Bumps',  'Blocks', 

%  'Doppler',  'Ramp',  'Cusp',  'Sing',  'HiSine', 

%  'LoSine',  'LinChirp',  'TwoChirp',  'QuadChirp', 

%  'MishMash',  ' WernerSorrows '  (Heisenberg), 

%  'Leopold'  (Kronecker) ,  'Piece-Regular'  (Piece-Wise 


Smooth) , 

%  ' Riemann ' , ' HypChirps ' , ' LinChirps ' ,  'Chirps',  'Gabor' 

%  ' sineoneoverx ' , 'Cusp2', ' SmoothCusp ' , 'Gaussian' 

%  'Piece-Polynomial'  (Piece-Wise  3rd  degree  polynomial) 

%  n  desired  signal  length 

%  Outputs 

%  sig  1-d  signal 

o. 

o 

%  References 

%  Various  articles  of  D.L.  Donoho  and  I.M.  Johnstone 

o. 

o 


if  nargin  >  1, 
t  =  ( 1 : n)  . /n; 
end 

if  strcmp (Name, ' HeaviSine ') , 
sig  =  4 . *sin  (4*pi  .  *t) ; 

sig  =  sig  -  sign(t  -  .3)  -  sign(.72  -  t); 
elseif  strcmp (Name ,' Bumps ') , 


pos  = 

[  -1 

.13  .15 

.23 

.25 

.40 

.44 

.  65 

.76 

.78 

•  81]  ; 

hgt  = 

[  4 

5  3 

4 

5  4  . 

.2  2. 

.  1  4  . 

.3 

3.1  5. 

.  1  4  . 

.2]  ; 

wth  = 

[  .005 

.005  . 

006 

.01  . 

.  01  . 

.03  . 

.01 

.01  . 

.005 

.008 

. 005]  ; 

sig  =  zeros  ( size  (t )) ; 
for  j  =1 : length (pos ) 

sig  =  sig  +  hgt (j)./(  1  +  abs ( ( t  - 
pos  (  j) )  . / wth ( j )  ) )  . A4; 
end 

elseif  strcmp (Name, ' Blocks ') , 

pos  =  [  .1  .13  .15  .23  .25  .40  .44  .65  .76  .78  .81]; 

hgt  =  [4  (-5)  3  (-4)  5  (-4.2)  2.1  4.3  (-3.1)  2.1  (- 

4.2)  ]  ; 

sig  =  zeros  ( size  (t )) ; 
for  j=l : length (pos ) 

sig  =  sig  +  (1  +  sign (t-pos ( j ))).* (hgt ( j ) /2 )  ; 

end 
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elseif  strcmp (Name, ' Doppler ') , 

sig  =  sqrt ( t . * ( 1-t )  )  . *sin  (  ( 2*pi*l . 05 )  ./(t+.05)); 

elseif  strcmp (Name, ' Ramp ') , 
sig  =  t  -  (t  >=  .37); 
elseif  strcmp (Name, ' Cusp ') , 

sig  =  sqrt (abs (t  -  .37)); 
elseif  strcmp (Name, ' Sing ') , 
k  =  floor (n  *  .31); 
sig  =  1  ./abs(t  -  (k+.5)/n); 
elseif  strcmp (Name, ' HiSine ') , 

sig  =  sin (  pi  *  (n  *  .6902)  .*  t); 

elseif  strcmp (Name, ' LoSine ') , 

sig  =  sin (  pi  *  (n  *  .3333)  .*  t); 

elseif  strcmp (Name, ' LinChirp ') , 

sig  =  sin  (pi  .*  t  .*  (  (n  .*  .500)  .*  t) ) ; 

elseif  strcmp (Name, ' TwoChirp ') , 

sig  =  sin  (pi  .*  t  .*  (n  .*  t ) )  +  sin((pi/3)  .*  t  .*  (n 

•  *  t )  )  ; 

elseif  strcmp (Name, ' QuadChirp ') , 

sig  =  sin (  (pi/3)  .*  t  .*  (n  .*  t.A2)); 

elseif  strcmp (Name, ' MishMash ') ,  %  QuadChirp  +  LinChirp  + 

HiSine 


sig  = 

sin 

(  (pi/3)  .* 

t  .  * 

( n  .  *  t . ' 

'2)  )  ; 

sig  = 

sig 

+  sin  (  pi 

*  (n 

*  .6902) 

.*  t)  ; 

sig  = 

sig 

+  sin  (pi  , 

.  *  t  . 

.*  (n  .* 

.  125  . *  t) ) ; 

elseif  strcmp (Name, ' WernerSorrows ' ) , 

sig  =  sin (  pi  .*  t  .*  (n/2  .*  t.A2))  ; 

sig  =  sig  +  sin (  pi  *  (n  *  .6902)  .*  t); 

sig  =  sig  +  sin  (pi  .*  t  .*  (n  .*  t)); 

pos  =  [  .1  .13  .15  .23  .25  .40  .44  .65  .76  .78  .81]; 

hgt  =  [  4  5  3  4  5  4. 22. 14. 3  3. 15. 14. 2]; 

wth  =  [.005  .005  .006  .01  .01  .03  .01  .01  .005  .008 

.005]  ; 

for  j  =1 : length (pos ) 

sig  =  sig  +  hgt (j)./(  1  +  abs ( ( t  - 
pos  (  j) )  . / wth ( j )  ) )  . A4 ; 
end 

elseif  strcmp (Name, ' Leopold ') , 

sig  =  (t  ==  floor (.37  *  n)/n);  %  Kronecker 

elseif  strcmp (Name, ' Riemann ') , 
sqn  =  round ( sqrt (n) ) ; 

sig  =  t  .*  0;  %  Riemann's  Non-dif f erentiable  Function 

sig ( ( 1 : sqn) . A2 )  =  1.  ./  (l:sqn); 

sig  =  real  ( if ft ( sig) ) ; 

elseif  strcmp (Name, ' HypChirps ') ,  %  Hyperbolic  Chirps  of 
Mallat ' s  book 

alpha =  15*n*pi/1024 ; 
beta  =  5*n*pi/1024; 

t  =  (1. 001:1: n+. 001) ./n; 

f 1  =  zeros  ( 1 , n) ; 

f 2  =  zeros  ( 1 , n) ; 

fl  =  sin (alpha ./(. 8-t )).*( 0 . l<t ).* ( t  <  0 . 68 ) ; 
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f 2  =  sin (beta. / ( . 8-t) )  . *  (0  .  l<t)  . *  (t<0 . 75) ; 

M  =  round ( 0 . 65*n) ; 

P  =  floor (M/ 4 ) ; 

enveloppe  =  ones (1,M) ;  %  the  rising  cutoff  function 
enveloppe ( 1 : P )  =  ( 1  +  sin (-pi/2+ (( 1 : P ) -ones ( 1 , P ))./ (P- 

1) *pi) ) / 2 ; 

enveloppe (M-P+l : M)  =  reverse (enveloppe ( 1 : P )) ; 
env  =  zeros ( 1 ,  n)  ; 

env (ceil (n/10 ) : M+ceil (n/10 ) -1 )  =  enveloppe ( 1 : M) ; 
sig  =  (fl+f2) ,*env; 

elseif  strcmp (Name, ' LinChirps ' ) ,  %  Linear  Chirps  of 
Mallat ' s  book 

b  =  100*n*pi/ 1024 ; 

a  =  250*n*pi/ 1024 ; 

t  =  ( 1 : n) . /n; 

A1  =  sqrt  (  (t-l/n)  . * ( 1-t ) ) ; 

sig  =Al.*(cos((a*(t)."2))  +cos((b*t+a*(t)."2))); 
elseif  strcmp (Name, ' Chirps ') ,  %  Mixture  of  Chirps  of 
Mallat ' s  book 

t  =  (l:n) ./n.*10.*pi; 

fl  =  cos (t . "2*n/1024 ) ; 

a  =  30*n/ 1024 ; 

t  =  ( 1 : n)  . /n . *pi; 

f2  =  cos  (a  .  *  (t .  "3 ) )  ; 

f2  =  reverse ( f  2 ) ; 

ix  =  (-n:n) ./n.*20; 

g  =  exp (-ix . "2*4*n/1024 ) ; 

11  =  (n/2  +  1 :  n/2+n)  ; 

12  =  (n/8+1 : n/8+n) ; 

j  =  ( 1 :  n )  /  n ; 

f3  =  g ( i 1 )  . *cos ( 50 . *pi . * j *n/1024 )  ; 
f4  =  g ( i2 )  . *cos  ( 350 . *pi . * j *n/ 1 024 )  ; 
sig  =  fl+f2+f3+f4; 

enveloppe  =  ones  (l,n) ;  %  the  rising  cutoff  function 
enveloppe ( 1 : n/ 8 )  =  ( 1  +  sin  (-pi/2+ (( 1 :  n/8 ) - 
ones  ( 1 , n/8 ) )  ./(n/8-l)*pi))/2; 

enveloppe  (7*n/8+l : n)  =  reverse (enveloppe ( 1 : n/8 )) ; 
sig  =  sig . *enveloppe ; 

elseif  strcmp (Name, 'Gabor'),  %  two  modulated  Gabor 
functions  in 


%  Mallat 's  book 


N  =  512; 
t  =  ( -N : N) *5/N; 
j  =  (1 :N) ./N; 
g  =  exp (-t . "2*20)  ; 

11  =  ( 2 *N/4+l : 2 *N/4+N) ; 

12  =  (N/4+1 :N/4+N) ; 

sigl  =  3*g ( i 1 ) . *exp ( i*N/l 6 . *pi . * j ) ; 
sig2  =  3*g ( i2 ) . *exp ( i*N/4 . *pi . * j ) ; 
sig  =  sigl+sig2; 

elseif  strcmp (Name ,' sineoneoverx ') ,  %  sin(l/x)  in  Mallat 's 


book 
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N  =  1024; 
i  =  ( — N+l : N) ; 
i  (N)  =  1/100; 

1=1./ (N— 1 ) ; 
sig  =  sin(1.5./(i)); 
sig  =  sig (513 : 1536) ; 
elseif  strcmp (Name, ' Cusp2 ' ) , 

N  =  64; 

1  =  ( 1 :  N)  .  /N; 

x  =  (l-sqrt(i))  +  i/2  -.5; 

M  =  8*N ; 

sig  =  zeros ( 1 , M) ; 
sig (M-l . 5 . *N+1 :M- . 5*N)  =  x; 
sig (M-2 . 5*N+2 :M-1 . 5 . *N+1)  =  reverse (x) ; 
sig ( 3*N+1 : 3*N  +  N)  =  . 5*ones ( 1 , N) ; 
elseif  strcmp (Name, ' SmoothCusp ' ) , 
sig  =  MakeSignal ( ' Cusp2 ' ) ; 

N  =  64; 

M  =  8*N ; 
t  =  (1  :  M)  /M; 
sigma  =  0.01; 
g  =  exp  (- . 5 . * (abs (t- 
. 5) . / sigma) . ^2) . / sigma . / sqrt (2*pi) ; 
g  =  fftshift (g) ; 
sig2  =  iconv (g ' , sig) ' /M; 
elseif  strcmp (Name, ' Piece-Regular ' ) , 
s igl=-l 5 *Make Signal ( ' Bumps ' , n) ; 
t  =  ( 1 : f ix (n/12 ) )  ./fix(n/12); 

sig2=-exp (4*t) ; 
t  =  ( 1 : f ix (n/7 )  )  ./fix(n/7); 

sig5=exp ( 4  *t ) -exp ( 4 ) ; 
t  =  (l:fix(n/3))  ./fix(n/3); 

sigma=6/ 40; 

sig6=-7  0*exp (-(  (t-1/2)  . *  (t-J/2 )  )/ ( 2*sigma~2 ) ) ; 
sig(l:fix(n/7) )=  sig6(l:fix(n/7) ) ; 

sig( ( f ix (n/7 ) +1 )  : f ix  (n/5 ) )=0.5*sig6( ( f ix (n/7 ) +1 )  : f ix  (n/5 ) ) ; 

sig( (fix(n/5)+l)  : f ix (n/3 ) )=sig6 (  ( f ix (n/5 ) +1 )  : f ix (n/3 ) ) ; 

sig( (fix(n/3)+l)  : f ix (n/2 ) )=sigl (  ( f ix (n/3 ) +1 )  : f ix (n/2 ) ) ; 
sig( (fix(n/2)+l) : (fix(n/2)+fix(n/12) ) )=sig2; 
sig( (fix (n/2) +2*fix (n/12) ) :  - 
1: (fix(n/2)+fix(n/12)+l))=sig2; 

sig (fix (n/2) +  2*  fix (n/12) +fix  (n/20) +1 :  (fix (n/2) +2*fix (n/12) +3* fix ( 
n/20 ) ))=... 

-ones (1,  fix (n/2) +2*fix (n/12) +3*fix (n/20) -fix (n/2) -2*fix  (n/12) - 
fix (n/20) ) *25; 

k=f ix (n/2) +2*fix (n/12) +3*fix (n/20) ; 
sig(  (k+1)  :  (k  +  fix(n/7) ) )=sig5; 
dif f=n-5*f ix (n/5 ) ; 

sig (5*fix (n/5) +1 :n) =sig (diff : -1 : 1)  ; 
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%  zero-mean 
bias=sum ( sig) /n; 
sig=bias-sig; 

el seif  strcmp (Name, ' Piece-Polynomial ' ) , 
t  =  (l:fix(n/5))  ./fix(n/5); 

sigl=2  0*  (t.A3+t.A2+4)  ; 
sig3=40* (2 . *t . A3+t)  +  100; 
sig2=10 . *t . A3  +  45; 
sig4=16*t.A2+8.*t+16; 
sig5=20* (t+4 ) ; 

sig6(l:fix(n/10)) =ones (l,fix(n/10) ) ; 
sig6=sig6*20 ; 
sig(l:fix(n/5) )=sigl; 

sig(2*fix(n/5) : -1 : ( f ix (n/5 ) +1) ) =sig2 ; 

sig( (2*fix(n/5)+l) :3*fix(n/5) )=sig3; 

sig(  (3*fix(n/5)+l)  :4*fix(n/5) )=sig4; 

sig (  (4*f ix (n/5) +1)  : 5* fix (n/5) ) =sig5 (fix (n/5)  :-l  :  1) ; 

dif f=n-5*f ix (n/5 ) ; 

sig ( 5*f ix (n/5) +1 :n) =sig (diff : -1 : 1)  ; 

%sig( (fix(n/20)+l) : (fix(n/20)+fix(n/10) ) )=- 
ones (1, f ix  (n/10 ) ) *20; 

sig((fix(n/20)+l) : (fix(n/20)+fix(n/10)))=ones(l,fix(n/10))* 

10; 

sig( (n-fix(n/10)+l) : (n+fix(n/20)- 
f ix (n/1 0 ) ) ) =ones  (l,fix(n/20) ) *1 50 ; 

%  zero-mean 
bias=sum (sig) /n ; 
sig=sig-bias ; 

elseif  strcmp (Name, ' Gaussian ') , 
sig=GWN (n,beta) ; 
g=zeros ( 1 , n) ; 
lim=alpha*n; 
mult=pi/ (2*alpha*n) ; 
g ( 1 : lim) = ( cos (mult* ( 1 : 1 im) ) ) . A2 ; 
g  (  (n/2  +  1 )  :n)=g(  (n/2)  : — 1 : 1 ) ; 
g  =  rnshift (g, n/2 ) ; 
g=g/ norm (g) ; 
sig=iconv (g, sig) ; 

else 

disp ( sprint f (' MakeSignal :  I  don*t  recognize 
<<%s>> ' , Name ) ) 

disp (' Allowable  Names  are:') 
disp ( ' HeaviSine ' ) , 
disp ( ' Bumps ' ) , 
disp ( ' Blocks ' ) , 
disp ( ' Doppler ' ) , 
disp ( ' Ramp ' ) , 
disp (' Cusp ') , 
disp ( ' Crease ' ) , 
disp (' Sing ') , 
disp ( ' HiSine ' ) , 
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o\°  o\°  o\° 


disp 

' LoSine ' ) , 

disp 

' LinChirp ' ) , 

disp 

' TwoChirp ' ) , 

disp 

' QuadChirp ' ) , 

disp 

' MishMash ' ) , 

disp 

' WernerSorrows ' ) , 

disp 

' Leopold ' ) , 

disp 

'Sing' ) , 

disp 

' HiSine ' ) , 

disp 

' LoSine ' ) , 

disp 

' LinChirp ' ) , 

disp 

' TwoChirp ' ) , 

disp 

' QuadChirp ' ) , 

disp 

' MishMash ' ) , 

disp 

'WernerSorrows ' ) , 

disp 

' Leopold ' ) , 

disp 

' Riemann ' ) , 

disp 

' HypChirps ' ) , 

disp 

' LinChirps ' ) , 

disp 

' Chirps  '  ) , 

disp 

' sineoneoverx ' ) , 

disp 

' Cusp2  '  ) , 

disp 

' SmoothCusp  '  ) , 

disp 

' Gabor ' ) , 

disp 

'Piece-Regular' ) ; 

disp 

' Piece -Polynomial ' ) ; 

disp 

' Gaussian ' ) ; 

end 


Originally  made  by  David  L.  Donoho. 
Function  has  been  enhanced. 


o. 

o 

%  Part  of  WaveLab  Version  802 
%  Built  Sunday,  October  3,  1999  8:52:27  AM 
%  This  is  Copyrighted  Material 
%  For  Copying  permissions  see  COPYING. m 
%  Comments?  e-mail  wavelab@stat.stanford.edu 

o. 

o 
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2, 


Kurtosistest 


kurtosistest 


implements  thebootstrap-based  hypothesis 
test  H:  kurtosis  of  the  input=3  against 
kurtosis  of  input~=3  and  performs  soft  and 
hard  thresholding  on  the  data 


SYNTAX  :  [ soft, hard, kur]  =  kurtosistest [data, kurlim] 

INPUT  :  data:  data  sequence  to  be  tested 

kurlim:  kurtosis  limit  not  to  be  exceeded 
thresholding 


OUTPUT 


soft :  soft  thresholded  data  s< 
hard:  hard  thresholded  data  s< 
kur:  kurtosis  value  used  in  soft  thresholding^ 


SUB  FUNC 


none 


Written  by  Hasan  E.  KAN 


'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S'S- 

oooooooooooooooooooooooooooooooooooooooooooooooo 


function  [varargout ] =kurtosistest (varargin) ; 

datain=varargin { 1 } ; 
kurlim=varargin { 2 } ; 

hardout=zeros (size (datain)  )  ; 

test=boottest (datain,  'kurtosis', 3,1,  .05,99,25); 
if  test==l 

sof tout=datain ; 
hardout=datain ; 
kur=-l ; 

else 

storage=[] ; temp=datain; 

kur=mean (bootstrp ( 99, 'kurtosis ' , datain) ) ; 
while  kur>kurlim 

if  length (storage ) > round (length (temp) / 3 ) 
kur=kurlim; 

else 

locations=f ind (abs (temp) ==max (abs (temp) ) )  '  ; 
storage= [ storage;  locations]; 
temp  (locations ) =0 ; 

kur=mean (bootstrp ( 9  9,  'kurtosis ' , temp) )  ; 

end 


end 
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sof tout=datain*abs (3-kur) / (kurlim-3) ; 
if  -isempty ( storage ) 

sof tout ( storage ) =datain ( storage ) ; 
hardout ( storage ) =datain ( storage ) ; 

end 

end 

varargout{l}  =  softout; 
varargout{2}  =  hardout; 
varargout{3}  =  kur; 
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APPENDIX  B.  SIMULATION  RESULTS 


Simulation  results  for  the  three  test  signal  types  are  presented  in  detail  in  this 
Appendix.  Descriptions  for  these  test  signals  were  given  in  Chapter  4.  Results  include 
performance  criteria  obtained  for  the  MSE  and  distance  measure  defined  in  Chapter  4  for 
SNR  levels  between  -6  to  6  dB,  for  one  hundred  trials  per  SNR  level.  Finally,  the 
original,  noisy  and  the  recovered  versions  of  one  trial  are  shown  for  SNR  values  -6  dB, 

0  dB  and  6  dB  to  provide  a  visual  representation  of  the  FFT-based  denoising  scheme. 

A.  FFT-BASED  DENOISING 
1.  Sinusoidal  Signal 

The  frequency  for  the  sinusoidal  signal  was  selected  randomly  for  each  of  the  100 
trials  so  as  to  generalize  the  simulation  results.  The  frequency  was  limited  to  avoid 
aliasing  or  DC  signals.  Minimum  and  maximum  allowed  frequencies  for  the  sinusoidal 
test  signal  in  the  simulations  were  2/5 1 2  Hz  and  241/512  Hz  respectively,  where  the 
sampling  frequency  was  selected  as  1  Hz.  Figure  13  shows  the  performance  curves  for 
FFT-based  denoising  scheme.  These  curves  indicate  that  the  thresholding  schemes 
perform  similarly  on  this  type  of  signal.  Figures  14  to  16  illustrate  noisy  and  recovered 
versions  of  a  sinusoidal  signal  at  SNR  values  of  -6  dB,  0  dB  and  6  dB  respectively. 
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Figure  13.  FFT-based  denoising  for  sinusoidal  signal  type. 


Original  signal 


Figure  14.  FFT-based  denoising  for  sinusoidal  signal  type,  SNR=-6  dB. 
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Original  signal 


Figure  15.  FFT-based  denoising  for  sinusoidal  signal  type,  SNR=()  dB. 
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Original  signal 


Figure  16.  FFT-based  denoising  for  sinusoidal  signal  type,  SNR=6  dB. 
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2.  Constant  Amplitude  Chirp  Signal 

Soft  and  hard  thresholding  perform  similarly  according  to  MSE.  However,  the 
distance  measure  shows  that  soft  thresholding  outperforms  hard  thresholding,  especially 
at  low  and  medium  SNRs.  This  result  can  be  seen  in  Figure  17.  Figures  18  to  20  illustrate 
noisy  and  recovered  versions  of  a  sinusoidal  signal  at  SNR  values  of  -6  dB,  0  dB  and 
6  dB  respectively. 


(a)  Distance  measure  for  constant  amplitude  chirp  signal 


(b)  MSE  for  constant  amplitude  chirp  signal 


SNR 

Figure  17.  FFT-based  denoising  for  constant  amplitude  chirp  type  signal. 
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Original  signal 


Figure  18.  FFT-based  denoising  for  constant  amplitude  chirp  type,  SNR=-6  dB 


Original  signal 


Figure  20.  FFT-based  denoising  for  constant  amplitude  chirp  type,  SNR=6  dB. 
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3.  Chirp  with  an  Amplitude  that  Increases  in  an  RC  Time  Constant 
Fashion 

The  results  shown  in  Figure  21  indicate  that  the  thresholding  schemes  have 
similar  MSE  performances  at  SNR  values  lower  than  -3  dB,  whereas  soft  thresholding 
performs  better  at  higher  SNRs.  According  to  the  distance  measure  d,  soft  thresholding  is 
better  at  all  SNR  levels.  Figures  22  to  24  illustrate  noisy  and  recovered  versions  of  a 
sinusoidal  signal  at  SNR  values  of  -6  dB,  0  dB  and  6  dB  respectively. 


(a)  Distance  measure  for  increasing  amplitude  chirp  signal 


SNR 

Figure  21.  FFT-based  denoising  for  increasing  amplitude  chirp  type  signal. 
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Original  signal 


Figure  22.  FFT-based  denoising  for  increasing  amplitude  chirp  type,  SNR=-6  dB. 
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Original  signal 


0  50  100  150  200  250  ,  300,  350  400  450  500 

Recovered  signal  with  hard  thresholding 


Figure  23.  FFT-based  denoising  for  increasing  amplitude  chirp  type,  SNR=0  dB. 


Original  signal 


0  50  100  150  200  250  ,  300  350  400  450  500 

Recovered  signal  with  hard  thresholding 

2 

0 

-2 


0  50  100  150  200  250  300  350  400  450  500 


Figure  24.  FFT-based  denoising  for  increasing  amplitude  chirp  type,  SNR=6  dB. 
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B.  WT-BASED  DENOISING 

Many  different  parameters  affect  the  performance  of  the  WT- based  denoising 
scheme  such  as  the  number  of  decomposition  levels,  processing  method  for  the 
approximation  coefficients  and  the  order  of  the  Daubechies  wavelet  used.  Results  show 
no  significant  difference  among  the  Daubechies  wavelet  orders  considered.  Only  best 
results  for  each  signal  are  illustrated  in  this  section.  Results  shown  include  MSE  and 
distance  measure  performances  for  SNR  levels  between  -6  to  6  dB,  where  one  hundred 
trials  are  used  for  each  SNR  level.  Noisy  and  recovered  versions  of  one  trial  are  shown 
for  SNR  values  -6  dB,  0  dB  and  6  dB  to  provide  a  visual  representation  of  the  WT-based 
denoising  scheme. 

1.  Sinusoidal  Signal 

Soft  thresholding  outperforms  hard  thresholding  at  low  and  medium  SNR  levels 
as  shown  in  Figure  25.  However,  it  does  not  perform  as  good  as  any  of  the  thresholding 
methods  used  in  the  FFT-based  denoising  case  for  this  signal  type.  Noisy  and  recovered 
versions  of  a  sinusoidal  signal  at  SNR  values  of  -6  dB,  0  dB  and  6  dB  are  shown  in 
Figures  26  to  28. 
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(a)  Distance  measure  for  sinusoidal  signal,  4  level  decomposition,  method  2 


SNR 

Figure  25.  WT-based  denoising  for  sinusoidal  signal  type,  4  level,  method  2,  with  Daubechies 

wavelet  orders  3,  4  and  5. 
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Figure  26.  WT-based  denoising  for  sinusoidal  signal  type,  4  level,  method  2,  Daubechies  wavelet  4, 

SNR=-6  dB. 
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14 


2. 


Constant  Amplitude  Chirp 


Results  in  Figure  29  show  that  soft  and  hard  thresholding  performances  are 
similar  for  this  signal  when  using  WT- based  denoising  scheme.  The  WT-based  scheme 
performs  better  than  the  FFT-based  scheme  on  this  signal.  Figures  30  to  32  illustrate 
noisy  and  recovered  versions  of  a  sinusoidal  signal  at  SNR  values  of  -6  dB,  0  dB  and 
6  dB  respectively. 

(a)  Distance  measure  for  constant  amplitude  chirp  signal,  4  level  decomposition,  method  1 


Figure  29.  WT-based  denoising  for  constant  amplitude  chirp  type,  4  level,  method  1,  with 

Daubechies  wavelet  orders  3,  4  and  5. 
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Figure  31.  WT-based  denoising  for  constant  amplitude  chirp  type,  4  level,  method  1,  Daubechies 

wavelet  4,  SNR=0  dB. 
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0  50  100  150  200  250  300  350  400  450  500 


Figure  32.  WT-based  denoising  for  constant  amplitude  chirp  type,  4  level,  method  1,  Daubechies 

wavelet  4,  SNR=6  dB. 
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3.  Chirp  with  an  Amplitude  that  Increases  in  an  RC  Time  Constant 
Fashion 

Soft  and  hard  thresholding  methods  perform  similarly  according  to  the 
performance  criteria  shown  in  Figure  33.  WT-based  denoising  scheme  outperforms  FFT- 
based  denoising  scheme  on  this  signal.  Figures  34  to  36  illustrate  noisy  and  recovered 
versions  of  a  sinusoidal  signal  at  SNR  values  of  -6  dB,  0  dB  and  6  dB  respectively. 

(a)  Distance  measure  for  increasing  amplitude  chirp  signal,  4  level  decomposition,  method  1 


-6  -4  -2  0  2  4  6 

SNR 


Figure  33.  WT-based  denoising  for  increasing  amplitude  chirp  type,  4  level,  method  1,  with 

Daubechies  wavelet  orders  3,  4  and  5. 
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Original  signal 


Figure  34.  WT-based  denoising  for  increasing  amplitude  chirp  type,  4  level,  method  1,  Daubechies 

wavelet  5,  SNR=-6  dB. 
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Figure  36.  WT-based  denoising  for  increasing  amplitude  chirp  type,  4  level,  method  1,  Daubechies 

wavelet  5,  SNR=6  dB. 
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APPENDIX  C.  COMPARISON  WITH  WAVELET  SHRINKAGE 

ALGORITHM 


The  proposed  denoising  schemes  are  compared  with  the  original  Wavelet 
shrinkage  algorithm  in  this  chapter.  Wavelet  shrinkage  algorithm  is  a  well-known 
denoising  scheme. 


A.  WAVELET  SHRINKAGE  ALGORITHM 

The  Wavelet  shrinkage  algorithm  was  introduced  by  Donoho  and  Johnstone  [16, 
17]  to  denoise  signals  embedded  in  additive  white  Gaussian  noise  with  unit  variance. 

This  algorithm  has  three  steps: 

a.  Wavelet  transformation. 


b.  Noisy  coefficients  suppression  by  applying  a  non-linear  thresholding 
technique. 

c.  Inverse  wavelet  transformation. 

1.  Choosing  a  Threshold  Value 

The  following  four  threshold  selection  schemes  are  available: 

a.  Stein’s  Unbiased  Risk  Estimator  (SURE)  Threshold 

The  SURE  threshold  value  is  derived  adaptively  for  each  decomposition 
level  by  minimizing  the  Stein’s  Unbiased  Estimate  of  risk  [18]  for  threshold  estimates. 

b.  Sqtwolog  Threshold 

This  method  uses  a  fixed  threshold  value  defined  as: 

T  =  J2\og(length(signal))  . 

c.  Heursure  Threshold 


This  method  is  a  mixture  of  the  preceding  two  methods.  It  uses  the 
Sqtwolog  threshold  at  low  SNR  levels  and  the  SURE  threshold  at  medium  to  high  SNR 
levels. 
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d.  Minimax  Threshold 

This  method  uses  a  fixed  threshold  that  gives  minimax  performance  for 
the  MSE.  The  minimax  principle  is  used  in  statistics  in  order  to  design  estimators,  which 
realize  the  minimum  of  the  maximum  mean  square  error  obtained  for  the  worst  function 
in  a  given  set.  So,  minimax  threshold  is  the  value  that  provides  the  minimum  MSE 
among  the  worst  threshold  values  for  the  detail  and  approximation  coefficient  sets. 

2.  Thresholding  Methods 

a.  Hard  Thresholding 

This  method  sets  the  coefficients  with  absolute  values  below  the  chosen 
threshold  to  zero. 

b.  Soft  Thresholding 

This  method  first  sets  the  coefficients  with  absolute  values  below  a  chosen 
threshold  to  zero.  Then,  it  shrinks  the  remaining  coefficients  using  the  relationship 

c  =  sign(c)[|c|—  7  ],  where  c  is  the  coefficient  to  be  thresholded,  c  is  the  thresholded 

coefficient,  and  T  is  the  chosen  threshold  value. 

B.  COMPARISONS 

Simulations  were  performed,  using  the  wavelet  shrinkage  algorithm  defined 
above,  on  the  same  signals  that  were  used  for  the  proposed  scheme.  However,  different 
SNR  values  were  obtained  by  using  constant  noise  power  and  different  values  for  the 
signal  power  in  this  experiment,  contrary  to  what  was  done  in  the  simulations  considered 
for  the  proposed  scheme  in  the  main  body  of  the  thesis.  The  Heursure  method  was  used 
to  compute  the  threshold  values  for  the  wavelet  shrinkage  algorithm.  A  new  simulation 
was  performed  with  the  proposed  scheme  under  the  same  SNR  conditions  as  well. 

Results  shown  in  Figures  37  to  39  indicate  that  the  proposed  WT-based  denoising  scheme 
performs  better  than  the  wavelet  shrinkage  algorithm  for  the  sinusoidal  test  signal. 

However,  the  wavelet  shrinkage  algorithm  outperforms  the  proposed  scheme  for  the  chirp 
type  test  signals.  It  should  be  noted  that  these  results  are  restricted  to  the  case  of  noise 
with  unit  variance. 
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1.  Sinusoidal  Signal 

Results  shown  in  Figure  37  indicate  that  the  MSE  performances  of  the  two 
denoising  algorithms  are  similar,  whereas  the  distance  measure  cl  shows  that  the  proposed 
denoising  scheme  outperforms  the  wavelet  shrinkage  algorithm  at  low  SNRs. 


(a)  Distance  measure  for  sinusoidal  signal,  3  level  decomposition 


SNR 

Figure  37.  Wavelet  shrinkage  algorithm  and  proposed  WT-based  denoising  scheme  performances 
for  the  sinusoidal  signal  type,  using  3- level  decomposition  and  Daubechies  wavelet  5. 
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2. 


Constant  Amplitude  Chirp 


Both  performance  criteria  show  that  the  wavelet  shrinkage  algorithm  outperforms 
the  proposed  denoising  scheme  for  this  signal.  This  can  be  seen  in  Figure  38. 


(a)  Distance  measure  for  constant  amplitude  chirp,  4  level  decomposition 


SNR 

Figure  38.  Wavelet  shrinkage  algorithm  and  proposed  WT-based  denoising  scheme  performances 
for  the  constant  amplitude  chirp  type,  using  4- level  decomposition  and  Daubechies  wavelet  4. 
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3.  Chirp  With  an  RC  Time  Constant-Like  Amplitude  Increase 

Both  performance  criteria  show  that  the  wavelet  shrinkage  algorithm  outperforms 
the  proposed  denoising  scheme  for  this  signal  as  shown  in  Figure  39. 


(a)  Distance  measure  for  increasing  amplitude  chirp,  4  level  decomposition 


SNR 

Figure  39.  Wavelet  shrinkage  algorithm  and  proposed  WT-based  denoising  scheme  performances 
for  the  increasing  amplitude  chirp  type,  using  4-level  decomposition  and  Daubechies  wavelet  5. 


It  should  be  noted  that  the  new  proposed  method  (based  on  kurtosis  and  Bootstrap 
method)  seem  to  be  inferior  to  the  methods  described  in  this  appendix.  But  the  methods 
in  this  appendix  need  additional  information:  the  variance  of  the  noise  is  unity,  which 
may  be  an  unreasonable  constraint  in  some  practical  applications. 
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